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SPECTRAL  ANALYSIS  OF  A UNIVARIATE  PROCESS  WITH 
BAD  DATA  POINTS,  VIA  MAXIMUM  ENTROPY  AND 
LINEAR  PREDICTIVE  TECHNIQUES 


1.  INTRODUCTION 

The  analysis  of  power  density  spectra  of  random  processes  via  maximum 
entropy,  linear  predictive,  and  autoregressive  techniques  has  attracted  much 
attention  recently,  especially  for  short  data  segments.  In  particular,  a good 
review  article  (reference  1)  recently  appeared  in  which  116  references  are 
listed  on  the  topic  of  linear  prediction.  Another  good  paper  on  this  method  of 
spectral  analysis  (including  a comparison  of  techniques)  is  available  in  refer- 
ence 2,  where  66  references  are  cited.  Additional  related  references,  that  this 
author  is  aware  of,  are  given  in  references  3 through  16  of  this  report.  The 
dose  links  that  exist  between  maximum  entropy  spectral  analysis  (MESA), 
autoregressive  spectral  analysis,  prediotive  error  filters,  noise-whitening  fil- 
ters, and  least-squares  model  building  are  pointed  out  very  well  in  reference 
14. 


The  purposes  of  this  report  are  to  review  and  interrelate  several  available 
techniques  for  spectral  analysis  under  different  states  of  knowledge,  for  equi- 
spaced  samples,  in  a consistent  notation;  collect  and  compare  the  techniques 
via  simulation  in  order  to  determine  the  best  available  technlque(s);  and  extend 
the  best  technique(s)  to  handle  the  case  of  bad  (or  missing)  data  points  and  com- 
pare them  via  simulation.  The  only  detailed  comparison  of  techniques  for  no 
missing  data  pointB  available  thus  far  in  the  literature  is  that  in  reference  2, 
where  the  Burg  technique  and  the  Yule-Walker  approach  are  compared.  Here 
we  will  extend  the  comparison  to  inolude  the  Burg  technique,  the  Yule- Walker 
approach,  an  unbiased  version  of  the  Yule-Walker  approaoh,  the  approximate 
maximum  likelihood  and  least-squares  approaches  of  reference  16,  the  auto- 
correlation and  covariance  approaches  of  reference  1,  and  an  extended  version 
of  the  covariance  approach.  (A  comparison  with  the  maximum  likelihood  tech- 
nique is  reserved  for  a future  report.)  Also,  we  will  compare  the  best  of  these 
approaches  for  the  case  of  bad  (or  missing)  data  points  and  present  FORTRAN 
programs  for  the  recommended  techniques. 

Throughout  this  report,  we  assume  we  are  dealing  with  equispaced  samples 
of  a stationary  zero-mean  random  process  x(t);  that  is,  x^sx(nA),  where  Ais  the 


1 
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sampling  interval  in  time.  In  section  2,  we  will  assume  that  the  correlation 
function  of  the  sampled  process,  jx^ . namely,* 


R,  s x x* 
k n n-k 


R-k’ 


(1) 


is  known  exactly  for  all  k,  and  shall  present  two  alternative  equations  to  deter- 
mine the  spectrum  of  |xnJ ; the  latter  of  the  two  equations  serves  as  a guide  to 
the  MESA,  linear  predictive,  and  autoregressive  approaches.  In  section  3,  it 
will  be  assumed  that  is  known  only  for  a limited  range  of  values  of  k,  and 
three  alternative  approaches  will  be  considered  and  shown  to  lead  to  identically 
the  same  spectral  approximation.  Next,  in  sections  4 and  5,  the  practioal 
problem  of  an  unknown  correlation  function  and  only  a finite  data  set  of  jxn|  , 
n « 1,  2,  . . . , N,  some  of  which  may.  be  bad,  will  be  addressed,  and  several 
candidate  techniques  for  spectral  estimation  will  be  presented.  Finally,  a com- 
parison of  the  techniques,  via  simulation,  will  be  conducted  and  conclusions 
drawn  regarding  the  best  available  technique,  both  with  and  without  bad  data 
points.  FORTRAN  programs  for  the  best  technique  for  both  situations  will  also 
be  presented. 


i 

I 


♦The  case  of  complex  samples  is  treated,  so  that  we  can  handle  complex 
envelope  or  complex  demodulated  processes.  Specialization  to  real  processes 
is  immediate , and  (1)  becomes  Rfc  = R-k.  An  overbar  indicates  an  ensemble 
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2.  CORRELATION  KNOWN  EXACTLY  FOR  ALL  ARGUMENT  VALUES 

Suppose  the  correlation  function  in  (1)  of  process  jxn|  is  known  for  all  k. 
The  standard  (double- sided)  definition  of  the  spectrum  of  |xn}  is  then  (see,  for 
example,  reference  14,  equation  (10)) 


Gx<f)  ■ 4 L Rk  exp(-i2nfkA),  If)  < (2) 

k=-CD 

Gx(f)  is  real  and  nonnegative,  but  need  not  be  even  in  frequency  f for  complex 
|Rk}- 

2. 1 LINEAR  PREDICTION  BASED  ON 
INFINITE  PAST 

Suppose  that  sample  values  xk_i,  xk_2»  • • • are  available  and  are  used  to 
linearly  predict  the  value  of  xk.  Then  the  one-step  predicted  value,  based  on  the 
infinite  past,  is  (for  a zero-mean  prooess) 


a x 
n 


k-n  ‘ 


(3) 


The  values  of  the  complex  predictive  filter  coefficients  j are  chosen  such 
that  the  one-step  prediction  error 


< 


k 


L 


n=  0 


a x 
n 


k-n 


(4) 


has  minimum  ensemble  average  magnitude-squared  value.  Figure  1 depicts  the 
Interrelationships . 


The  ensemble  average  magnitude- squared  error  is,  employing  (1),  given  by 

CD 


£c  U = T.  a a*  R 

Ik  ^.mnn-m 

m,n=  0 


(5) 


For  a minimum,  we  first  compute  (see  reference  17,  appendix  A) 
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Figure  1.  Block  Diagram  of  Predictive  and  Whitening  Operations 


■s..  £ *,  , ,jii 

* n l-m  m 


(8) 


da  f m=*  0 


and  set  It  equal  to  zero,  obtaining  the  optimum  predictive  filter  coefficients 
aB  the  solution  of  the  set  of  equations** 

1C  R/  m ^ 1 C2  =a  =-1) . 

JT-n  »-  m m o o 

m=  u 

The  minimum-error  sequence  jTjJ  then  possesses  correlation 


E.  S 7.7*  m £ * & X.  X* 

J k k"J  m,n=  0 m n k"m  k-^_n 


CD 


£ s a*  r 

• m n 

m,n=  0 


CD 


£ 

m"n  "J  +n-m  J^n  "n  ‘’J+n-m  m ' 


E *n*  E 8, 

n=  0 m=  0 


(7) 


(8) 


’•‘The  same  result,  (7),  can  be  obtained  by  setting  the  partial  derivatives  of 
E,  with  respect  to  the  real  and  imaginary  parts  of  a^,  equal  to  zero. 
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where  we  have  employed  (4)  and  (1) . Now  the  Innermost  sum  on  m In  (8)  is  0 
for  J + n £1,  by  (7).  And  if  j £ 1 in  (8),  then  j + n £ 1 sinoe  n 2 0 in  the  outer- 
most sum  in  (8).  Therefore,  Ej  = 0 for  j £1.  Also  sinoe  E_j  = EJ\  we  have 

Ej  - 0 for  j f 0 ; (9; 

that  is,  the  minimum-error  sequence  {7^}  is  uncorrelated  and  therefore  pos- 
sesses a white  speotrum.  The  linear  filter  characterized  by  coefficients 
{ } 0 is  a whitening  filter;  see  figure  1. 

The  correlation  of  { for  zero  time  delay  is  the  power  of  the  minimum 
error  and  is  given  by 


19  UJ  00 

- £ & £ R ** 
n-0  n m«0  n"m  1 


£ R * “ R - E ***  * , 

° m = 0 ~m  m 0 m-1  m m 


where  we  have  used  (8),  (7),  and  (1).  The  spectrum  of  17.  \ is  therefore 
(using  (9))  Kl 


00 

G-(f)  ®A  £ E e:cp(-i2nfjA)  ■ AE  , If  I < ~ , (11) 

j=  -0D  •*  O 2A 

which  is  white,  as  mentioned  above. 

But  sinoe  the  error  sequenoe  is  given  by  a linear  transformation  of  process 
l xk(  according  to  (4)  and  figure  1,  the  spectrum  of  {7^1  is  given  by  the  stand- 
ard linear  filter  relation 


where 


G?(f)  = |A(f)|2G^(f)  , 


A(f)  = £ S exp(-i2nfnA) 
n=  0 
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Is  the  transfer  function  of  the  whitening  filter  and  is  assumed  to  be  b table.* 
Combining  (11)- (13),  we  obtain  an  alternative  expression  to  (2)  for  the  spectrum 
of  jxk!  as 


G (f)  « 
x' 


4E 
0 


co 

2 a exp(-i2*rfn4) 
na  0 


(14) 


Given  the  correlation  values  {R^} , utilization  of  (14)  requires  solution  of  the 
set  of  equations  in  (7)  for  the  filter  coefficients  {S^}  and  subsequent  substitu- 
tion in  (10)  and  (14).  Although  this  is  not  a practical  alternative  to  (2)  in  this 
case,  it  does  serve  to  indicate  that  there  is  posBlbly  some  potential  in  the  idea 
of  determining  predictive  filter  coefficients  to  minimize  the  average- magnitude- 
squared  one-step  prediction  error  and  thereby  obtain  a white  speotrum;  this 
idea  will  prove  to  be  quite  fruitful  later  on. 

As  an  aside,  if  we  allow  aM^  4 0 in  (3)  and  minimize  |«h)2i  we  find  4 0, 
although  Ej  » 0 for  j >2.  Thus,  the  minimum-error  sequence  would  not  be 
white,  and  a convenient  expression  like  (14)  would  not  result. 


It  should  also  be  noted  that  the  orosscorrelations  between  the  minimum- 
error  sequence  |«jj  and  all  past  values  of  the  input,  |xjt[ , are  zero;  this 
follows  by  use  of  (4),  (1),  and  (7). 


2.  2 LINEAR  PREDICTION  BASED  ON 
INFINITE  FUTURE 


If  sample  values  xj(+1,  xk+2»  • • • 0X6  available  and  are  used  to  linearly 
"predict"  the  value  of  x^  according  to  a backward  regression  (that  is, 
combine  future  values), 


a*  x 
n 


k+n' 


(15) 


x _n 

*That  is,  _S  n tf  z has  all  its  poles  inside  the  unit  circle,  O,  in  the  com- 
- . n« u n 

plex  z -plane. 


• * ... 


; A 


t’J-V,  ' «*■»=«  • • * . r • - • -**  ■ * - - ' * « 
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then  the  one- step  error 


•k,*k-xk-  s AY, 

n*  u 


has  average  magnitude- squared  value 


m,n«0 


a a*  R , 
m n n-m 1 


whloh  is  identloal  to  (5).  Thus,  the  same  optimum  filter  coefficients  in  (7)  that 
minimized  (5)  would  also  minimize  (17).  The  minimum-error  sequence  in  (16) 
would  also  be  white,  and  an  expression  for  the  speotrum  of  jx^}  identical  to 
(14)  would  result.  The  point  of  this  result  is  that  an  equivalent  expression  for 
the  speotrum  of  |x^|  is  obtained  by  the  backward  regression  (15),  rather  than 
the  forward  regression  (3)  of  the  preceding  subsection.  This  idea  will  prove 
usefhl  later  when  we  have  to  deal  with  finite  data  sets  and  unknown  correlation 
functions. 

I, 

The  orossoor relations  between  the  minimum- error  sequenoe  and  all  future 
values  of  the  input  are  zero;  this  follows  by  use  of  (16),  (1),  and  (7). 


2.3  LINEAR  INTERPOLATION  BASED 
ON  INFINITE  PAST  AND  FUTURE 

If  we  attempt  to  combine  the  approaches  of  the  previous  two  subsections, 
we  are  led  into  considering  linear  interpolation  according  to 


x.  = £ ax, 
k n k-n 

n»-oo 

n/0 


The  error 


L Vk-„ 

n=  -oo 


<vl) 


, •*  - .Ik  t »* 
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has  average  magnitude-squared  value  > . , , u 

" »«  Kf  ” vs  v»*  ■ <20) 

m,  n**  itn 

using  (1).  Setting  AE/da*  ■ 0 for  0,  we  obtain. for { the  optimum  filter  qoeffli 
oients 


' ! "i. 


“ i-m  m 


'Q-ao--l).  4/  (21) 


m®  -<u 

lie-  ■ • ■••■■i,. 


There  follows,  by  use  of ' (1),m  . ..  ■ > , , . . ,,  ,, 

l*  1 ••  •<•  ill  ! j.  ! ■ ' • 1 "I  •.  .1  . , .... 

■ ' ■ . ' ’ll  ■ "!•  I VI.  I V . . .11  ' 

'.I  "•  ' . v ■ i ''  . : • r,'  hi  ■ 


1 1 * ; i i'  - i ” 1 1 , i s i i . i 

. !,•••!  •! . . . 


i 'ii’4  . » ■ , i 1 1-  , i ' ■ t ..  ■■  , 

■ i ; : n » . * i . 1 ■ i i » ' • i [ [ i ' 1 1 \ . .. . j . 

1 1 iI'H.* y} 

.i  -I,...  (22) 

’•■111  ■ . > ,,  f I..  .11 


r Tt)e  correlation  of  the  minimum+jarror  sequei^oe  is  novlr  1 v". 


E4  a ST  7*  ~ =»  T If  3*  Rix 
J k k-j  m n j + n-m 

. ■ ■ it-  m, n» -«>  r,  . . •• 


i i : i f 


0)0)  O) 

« T a*  y K R.  , * V £ It  R =>  - * E , (23) 

n +-<  m j+n-m  . +j  i A*r.  * an'^ai".!  . J o’  V 

n=>-«>  m = -w  * * m=  -u> 

• I . • ! : ■ i ; * 1 1 ■ 

where  we  have  employed  (19),  (1),  (21),  and  (22),.  It  is  generally  nonzero  for 
j ^ 0^  *The  spectrum  of1  the  minimum-error  sequence  is 'therefore 

U>  1 

Gy(f)  - - AEo  £ a exp(-i2irfjA)  = - AEoA(f),  If  I < ~,  (24) 

j=-0U  J 

• I 

where  we  have  used  (23)  and  assumed  A(f)  to  be  Btable.  This  Qpectrum  is  not 
white;  in  fact,  employing  (12),  (24)  can  be  expressed  as 

A2E2 

Gf<»  - gT®  ; 111  < S-  <26> 

x' 

which  is  the  inverse  of  the  input  Bpeotrum . 
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If  we  instead, eliminate  G y(f)  from  (ia>  and  (24),  we  obtain  an  expression, 
for  the  input  spectrum  terms  of  filter  A(f)  (in  (24)  as  , 

■i  -;1 if  ae  . . ' ' ■ 

, ■ 1 1 1 1 ■ i ! " A*(f):‘''A(f)  • 2A:  • 


the  realness  of  A(f)  follows  from  (22). 


There  is  an  uncorrelated  property  between  the  minimum  orror  and  the  input 
in  the  present  oase  also.  Namely,  the  orossoorrelation  between  the  minimum- 
error  sequence  and  the  input  is 


V VI. 


CD 


E X 


n=  -a) 


n k-n  k-J 


t 


Am  -OD 


Z H. 
n J-n 


E 5 1 1 
o "oj 


(27) 


using  (19),  (1),  (21),  and  (23).  ThuB,  the  minimum-error  sequence  is  uncorre- 
lated with  all  past  and  future  values  of  the  input  exoept  at  the  same  time  instant. 
The  cross-sp#otinim  is 

OD  , 

Gj-^f)  5 * £ Cj  exP(-i2rrfjA)  - *Eo  , If  I < (28) 

j =5  — OJ  * J 


which  is  white. 

Although  (26)  and  (21)  offer  an  alternative  to  (14)  and  (7)  in  the  present  oase 
of  known  correlation  funotion  jlt^  j , it  suffers  in  the  practical  oaBe  of  unknown 
correlation  and  a finite  data  set,  by  virtue  of  the  estimate  of  the  real  denomina- 
tor of  (26)  going  through  zero  (or  being  complex  if  (22)  is  ignored)  at  some 
values  of  f . This  is  not  a significant  problem  for  (14)  since  both  the  real  and 
imaginary  parts  of  the  estimate  of  (13)  must  simultaneously  equal  zero  there, 
in  order  to  constitute  a problem. 

Another  important  practical  drawbaok  of  this  interpolation  approach  is  that 
ensemble  average  |«]J*  would  probably  be  approximated  by  | K!2* where 

the  sum  is  conducted  over  those  values  of  k at  which  a meaningful  value  of 
error  can  be  formed  for  a segment  of  a single  member  function  of  an  en- 
semble. But  since  the  minimum-error  sequence  is  not  unoor related  in 

this  case  (see  (23) ),  minimization  of  ^ |*kl^  *or  a sin#e  member  funotion 

segment  is  not  synonymous  with  minimization  of  |«k|2  ’ rather,  the  minimiza- 
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tion  of  £ jf^j2  will  spuriously  Involve  correlation  between  adjacent  terms 

which  are  not  inoluded  in  ) ejj.) 2 and  whioh  will  bias  the  filter  ooeffioienta . 
Several  simulation  runs  (on  real  data)  confirmed  this  oonoluslon  by  yielding 
severely  biased  (and  negative)  estimates  of  speotrum  Gx(f) , even  when  (22)  was 
taken  Into  aooount.  ^Accordingly,  the  interpolation  approaoh  was  dropped  from 
further  consideration. 


I 


10 
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3.  CORRELATION  KNOWN  EXACTLY  FOR  A LIMITED 
RANGE  OF  ARGUMENT  VALUES 

In  this  section,  Rk  of  (1)  is  assumed  to  be  known  exactly  for  Ik!  < p and 
unknown  for  Ikl  > p,  Since  we  are  unable  to  compute  the  exaot  spectrum  Gx(f), 
given  by  (2),  in  this  oase,  a different  approach  involving  approximation  Gx(f) 
is  required.  Three  different  techniques  will  be  considered  and  shown  to  yield 
identically  the  same  approximation  to  Gx(f). 

3.1  MAXIMUM  ENTROPY  SPECTRAL 
ANALYSIS  (MESA) 

The  method  in  this  subsection  was  originally  given  in  reference  18  and 
elaborated  upon  in  reference  19.  We  begin  with  (2)  and  note  that 

-L 

2A 


l df  Gx(f)  cxp(12fffkA)  ■ I df  Gx (f)  exp(i2nfkA)  • R^ 

\ * K/L 


(29) 


1 

'21 


We  wish  to  approximate  Gx(f)  by  a real  nonnegative  function  G(f)  such  that  its 
entropy  (reference  18,  equation  (1) ) 


l 


A | dfin  G(f) 
'1/A 


(30) 


1b  maximized,  subject  to  the  integral  oonBtraintB 

df  G(f)  exp(i27rfkA)  » Rk,  Ikl  5 p 


J 


1/A 

To  this  aim,  we  form  the  quantity 


Q s f dfin  G(f)  - £ f 
i/A  k=-P  Jl/A 


(31) 


df  G(f)  exp(i2trfkA) , 


(32) 
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where  Lagrange  multipliers  p_ka  , because  of  the  restriction  R_k  » Rk , as 
shown  in  (1).  We  perform  a variation  of  (32)  aooordlng  to 


Q+IQ 


df  jin 


[Go(f)  + «e(f)1  - E v(  df[Q0<f)+*^f)] 

■ >■  ■ ■ k«*-p  '""{fo  ■ ■ ■■■'■ 1 1)1 


. ■ l , ' , I : 


exp(i2»fk0) , 
(33) 


where  G0(f)  is  the  ’'optimum"  approximation  to  Gx(f)  undeir  criterion  (30),  and'  ' 
obtain,  upon  setting  ■ , , 


c*(Q  + 4Q) 


■ 0 at  « » 0 , 


the  relation 


G (f)  

ow  p 


■■  w<t' 


X nk  exp(i2»rfkA) 
k«-p 


GQ(f)  is  real  since  M_k  » nk . Since  it  is  also  to  be  nonnegative,  we  can  express 


where 


7(f)  ■ £ ak  exp(i2nfka),  Ifl  < jji 
k»  0 


and  whore  y(i)  has  no  zeros  in  the  upper-half  complex  f-pl^ne;  that  is,  poly- 
nomial 


B(z)  a £ “t  z 
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A 

l 


has  no  zeros  InBlde  the  unit  bl^elcS,  O,'  in  the  ddinplek  z~plane.  A prbof  that  ri 
B(z)  in  (38)  has  no  zeros  inside  O is  given  in  reference  11,  page  7,  for  exam- 
ple. * Specifically,  it  is  shpwn  that  B(l/z)  has  all  its  poles  and  zeros  inside  O; 
that  is,  B(l/z)  is  minimum  phase.  , , ; , 


In  order  .ip  determine  the,, constants  in  (37).  ,we  1expre^s,(3,6)  as, . 


. i nit!  n 

• * i.) 

• j ‘ " i 


•■'i » n 


I I 

Mill..  : i 


(39) 


(We  could  equally  well  have  multiplied  by  7(f)  in  the  following.)  Therefore,  for 
all  values  of  jJ , 1 


1 


df  Go(f)7*(f)  OXp(i2rrdA) 


/A 


f 


1 


df  -j.  exp(12*f^)..  (*0) 


1/A 


But  using  (37),  this  dan  be  expressed  as 
P 

I 

k<°  0 


£ a*  j"  df  Go(f)  exp(i2nf(j?-k)A)  * ( ' df  , all.fi. 

raft  » j.  *5  a «-**■*  


(41) 


1/A 


i/a  , £ «k  bxpd^fffKh)  ,, 
k=0 


Now  if  jfi  is  an  integer  in  the  range  [0,  p] , tho  integral  on  the  left  side  of  (41)  is 
equal  to  R|>k  (via  (29) ) fo^  any  value  of  k in  its  range  [0,  p] ; this  is  where  the 
constraints  are  employed.!  Therefore,  we  .have  fot*  integer  i, 


E R»-k 

kHO 


(42) 


where 


i ,i 


b' a J 


df. 


exp(i2nflA) 


, O' < i < p '." 


(43) 


1/A  ^ “k  exp(i2trfkA) 

k=0 


*The  proof  is  couched  in  terms  of  the  recursive  solution  of  (46)  presented 
in  appendix  A. 


T'3 


.......... 


.HIV  • 


m 
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In  (43),  letting  z « exp (12*14)  and  using  (38),  we  have 


i}.  & £ QOl  JL-  o < j2  < p 

* 12*4  7 z B(Z)  1 55  * P * 


where  $ denotes  counterolookwise  integration  around  the  unit  olrole  O in  the 
complex  z-plane,  Now  B(z)  has  no  zeros  inside  O by  oonBtruotion.  Further- 
more, B(z)  can  have  no  zeros  on  0,  for  then  y(f)  would  be  zero  for  some  f,  and 
G0(f)  would  possess  infinite  power,  oontradloting  R0<  Then  (44)  yields 


and  (42)  becomes 


k=>  0 

This  is  p+1  linear  equations  in  p+1  unknowns . * 
Now  let  correlation  matrix  R be  defined  as 


Ro  R-1 


R,  R 
1 o 


and  define  two  column  matrices 


I = tl  0 0 ...  0]T,  •=  [a  a . ...  a 1T 

0 1 D 1 


♦The  recursive  solution  of  (46)  is  presented  in  appendix  A. 
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i 

i 

t 


B is  Hermitian,  Toeplitz,  and  nonnegative  definite.  Then  (46)  oan  be  expressed 
as 


R «* 


Aa , 


(49) 


with  solution 


«*  = - 1 r"1  > 


(60) 


Now  let  the  inverse  matrix 


R-1  => 


c o , 
oo  ol 


cn 

lo  11 


po 


op 


pp 


(61) 


Then  (60)  and  (48)  yield 

.1/2 

a*  - — — • c ,[a|2--ro  , a - (——)  exp(ifl)  , (52) 

o oo  I o'  A oo  o \ A / ' ' 

where  0 Ib  an  arbitrary  real  constant.  (c00  is  always  real.)  Utilizing  this  re- 
sult and  (15)  in  (50),  there  follows 


kn 

a*  a — - — ■ - axnf-  i n\ . f)  < W <■  n . 
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I 


and  (37)  becomes 


■y(f) 


^7=T  E cko  exP(i2fffkA)«  lf  1 1 < h ' 
ka0 


(54) 


Finally,  using  (36),  the  "optimum"  speotrum  (called  the  maximum  entropy 
speotrum)  is 


Ac 


G0(f)  - 


oo 


P 

I 

k=  0 


£ Oko  exp(-  12rrfk  A) 


(55) 


Equation  (55)  gives  the  maximum  entropy  spectrum  in  terms  of  the  first 
column  of  the  inverse  of  the  correlation  matrix  R of  available  known  correlation 
values;  see  (47).  The  forms  of  (56)  and  (46)  are  similar  to  those  encountered 
earlier  in  (14)  and  (7),  respectively;  see  also  appendix  A,  The  maximum  value 
of  the  entropy  defined  in  (30)  is  evaluated  in  appendix  B and  is  given  by  Jin  (a/cqo). 


Substitution  of  (53)  in  (38)  yields 


B(z) 


exp (-ie) 


P 


E 

k=0 


(56) 


P 


Thus,  investigation  of  the  zeros  of  B(z)  depends  on  the  polynomial  2 c^  z^; 

k=  0 

it  must  have  no  zeros  inside  the  unit  circle  0.  But  if  we  combine  (46)  and  (53), 
we  can  write  that 


P cko 

E R,.kr  = -R/’  ^ 

k=  0 oo 


(57) 


Now  reference  1,  page  567,  declares  that  all  the  zeros  of  2 c^n  z“K  must  lie 

k=  0 

inside  O since  R is  a correlation  matrix.  Therefore,  polynomial  B(z)  has  no 
zeros  inside  0. 
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3.2  LINEAR  PREDICTIVE  FILTERING 

Here,  as  in  the  previous  subsection,  the  available  information  is  knowledge 
of  Rfc  for  Iki  < p.  A linear  one-step  prediotlon  of  xfc,  based  on  the  past  p values, 
xfc-i,  . . . , x^_p,  is  to  be  aooompllBhed  with  minimum  average  magnitude- 
squared  error;  see  figure  1.  Now,  however,  instead  of  (3),  we  have  for  the 
predicted  value  the  finite  sum* 


P 


£ 

n=  1 


a x 
n 


k-n  * 


The  instantaneous  error  is 


. a 
2 X, 


- X, 


P 

£ 

n°0 


a x. 
n k-n 


(ao  = -1) 


(58) 


(59) 


(Equations  (58)  and  (59)  constitute  stable  digital  filters  regardless  of  the  choioe 
of  coefficients.)  The  ensemble  average  magnitude-squared  error  is 


E 5 


P 

£ 


m,n=  0 


a a*  R 
m n n-m 


• HR  ■ 


where  we  have  used  (1)  and  (47)  and  defined 


(60) 


« = [a  a...,  a)1.  (61) 

o 1 p' 

We  now  wish  to  minimize  E by  choioe  of  filter  coefficients  . How- 
ever, we  have  the  constraint  on  a0  in  (59);  this  can  be  expressed  mathematically 
as 


J «H  • = - 1 , (62) 

where  I is  defined  in  (48).  In  order  to  minimize  (60)  subject  to  (62),  we  form 
the  quantity 


* P * 

The  more  general  form  including  ^ x^_n 


is  not  considered  here. 
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H _ . H .*  T_* 

- (■  - XR"1*)11  R(s  - XR-1*)  - IXI2»Hr“1«,  (03) 

where  R"1  is  defined  in  (51).  Since  R is  nonnegative  definite,  being  a correla- 
tion matrix,  (63)  is  minimized  by  the  choice  of  ooeffloients 

i . 

• ■ X R 1 * . (04) 

l 

j The  Lagrange  multiplier  X is  obtained  by  substituting  (64)  in  constraint  (62),  and 

| using  (51)  and  (48): 

i 

I 

! .*  -H  -1  . 1 

XI  R i a - 1,  » = - . (06) 

°oo 

Then  (64)  yields 

Cko 

\ “ o ’ 0 ~ k - p * (<*«) 

oo 

The  minimum  value  of  the  error  power  is  found  by  utilizing  (64)  and  (86)  in 

(60): 

E =|7.|2=  |XI2lHR"lRR"1|  - I XI2  o =>  -i-,  (67) 

o I k I oo  o ’ ' ' 

oo 

where  is  the  minimum-error  sequence  obtained  by  employing  (66)  in  (59). 

(A  re  ivirsion  relation  for  E^)  is  presented  in  (A-7);  it  can  be  started  with 
i l/°£-*o.>  Notice  from  (67)  that  cqq  must  be  positive,  for  non-negative 

definite  R . 

The  transfer  function  of  the  optimum  error  filter  from  input  x to  output  7 
in  figure  1 is,  from  (59)  and  (66), 
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A(f)  ■ 2 It  exp(-i27rfkA) 
k-0 


1 P i 

- — L °ko  exP(_l2nfkA)»  lfl  < 
oo  ka  0 


2A 


(68) 


Furthermore,  the  speotra  in  figure  1 are  related  by 


G —(f)  = i A(f)|  * 0 (f)  . 


(69) 


Now  let  us  assume  that  the  spectrum  of  the  minimum -error  sequence  is 
white  over  the  band  (-;rjr,  rr);  this  is  in  line  with  the  property  (11)  which  held 

UUk  £ A 

for  the  case  when  the  infinite  past  was  available.  Then  we  say 


lfl  < -- 
2A 


(70) 


where  we  have  used  (67) . Substitution  of  (68)  and  (70)  in  (69)  yields  the  linear 
predictive  spectrum  approximation  to  the  input  Bpeotrum  aooording  to  the  defi- 
nition 


G (f)  a 
x'  ’ 


Gy(f) 


a 


p 

53  0.  exp(-i2rrfkA) 
k»  0 


(71) 


This  is  identical  to  the  approximation  (55)  obtained  by  MESA.  It  is  critioally 
dependent  on  the  assumption  that  the  speotrum  of  the  minimum-error  7 in  fig- 
ure 1 io  white. 

Since  (71)  is  identical  with  the  maximum  entropy  spectrum,  (56),  it  must 
follow  that 


df  Gx(f)  exp(i2irfkA) 


R^  for  Ikl  < p j 


(72) 
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that  is,  although  not  specified  In  the  current  approaoh,  the  correlation  funotion 
formed  from  the  linear  predictive  spectrum  6rx(f)  in  (71)  has  the  same  values  at  \ 
kd  for  Ik  I <p  as  the  known  correlation  values  {Rjj}  . 

The  implications  of  the  assumption  (70)  of  a white  spectrum  for  the  mini- 
mum error  are  investigated  in  appendix  C.  It  is  shown  that  the  crossoorrela- 
tlon  funotion  between  input  x and  output  t of  figure  1, 

°<  5 V£-v 


must  then  satisfy 


0/ 


l/onn 

00 

I 

o ,i  >1, 


(74) 


that  is,  minimum-error  sequenoe  is  assumed  uncorrelated  with  all  the 

past  values  of  the  Input.  It  is  also  shown  that  the  unknown  correlation  values 
Rk  for  k > p can  be  approximated  according  to 


R, 


P o 

no 
. o 

n»l  oo 


R 


k-n 


P 

£ 

n=l 


VW 


k > p + 1 . 


(75) 


This  recursion  relation,  starting  with  known  values  Rj_,  . . . , Rp,  can  be  con- 
sidered to  be  an  extrapolation  of  the  known  correlation  values  into  regions 
where  they  are  unknown.  Equation  (75)  is  Bhown  in  appendix  D to  be  a stable 
recursion  when  B(z)  of  (56)  has  no  zeros  inside  O;  this  property  has  been  dis- 
cussed under  (38),  (56),  and  (A-9).  It  can  also  be  shown  that  Fourier  transfor- 
mation of  the  extrapolated  correlation  approximates  yields  precisely  (71).  It  is 
interesting  to  note  that  (75)  lias  the  same  form  as  the  prediotive  equation  (58)  for 
individual  data  values. 


Since 


df  6 (f)  exp(i2?rfkA) 


(76) 


is  the  autocorrelation  at  delay  kh,  it  is  given  by  (72)  for  Ikl  <p,  and  by  (75)  for 
k > p + 1,  where  the  latter  correlations  arc  extrapolated  values.  This  follows 


ft 
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from  setting  67(f)  white  and  choosing  Gx(f)  by  (71),  aooording  to  the  analysis  In 
appendixes  C and  D. 


If  sample  values  xk+i*  . . . , were  used  to  linearly  "predict"  x^  ac- 
cording to  backward  regression 


n«  1 

the  one-step  error  «k  = xk  ~ xk  has  average  magnitude- squared  value 

P 


(7  7 A) 


E=  l«J  - £ a <Rn 

1 1 ' m n n- 

m,n=  0 


m 


(ao=-D, 


(77B) 


| which  is  identical  to  (80).  Thus,  the  same  optimum  filter  ooefficients  in  (66) 

| that  minimized  (60)  would  also  minimize  (77B),  and  an  approach  similar  to  that 

( above  would  yield  a spectral  approximation  identical  to  (71).  The  equivalence 

of  the  results  of  this  backward  regression  to  that  of  the  forward  regression  in 
(58)  will  prove  useful  later  when  we  deal  with  finite  data  sets  and  unknown 
correlation  functions. 

3.3  ALL-POLE  DIGITAL  FILTER 
MODEL 

The  available  information  about  process  {x^}  is  the  same  as  In  the  previ- 
ous two  subsections,  namely,  knowledge  of  R^  for  lkl<  p.  Consider  a sampled 
autoregressive  process  jy^j  in  steady  state  generated  by  a stable  all-pole 
digital  filter,  H(z),  excited  by  discrete  white  noise  {w^}  ; see  figure  2,  The 
noise  is  characterized  by  correlation 


Wi_  w*  = b , all  n , (78A) 

k k-n  no 

with  no  loss  of  generality,  and  has 

G (f) 
w 


spectrum 

=>  A , Ifl  < ~ . (78B) 
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Figure  2.  Generation  of  Ail-Pole  Process 


The  digital  filter  is  characterized  by  a p-th  order  autoregressive  relationship, 


P 

yk-n  Wk  * 

n»0 

(79) 

with  transfer  function 

H(Z)  » 

£ V 

n=  0 

_ f k n 

(80) 

We  are  going  to  choose  digital  filter  coefficients  |j3n|  so  that  autoregres- 
sive process  {y^}  has  the  same  correlation  values  as  process  jxjc}  , up 
through  order  p;  that  is,  we  will  set 


y.  y*  = R for  Ini  < p. 
k k-n  n - * 

Then  the  spectrum  of  autoregressive  process  jy^  , given  by 


G (f)  = G (f) 
y'  w 


H 


(exp(i27rfA))j 


£ Pn  exp(-  i2nfnA) 
n=  0 


?m<^' 


(61) 


(82) 
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will  be  used  aa  an  approximation  to  the  speotrum  of  {xfc}  • The  epeotral  rela- 
tion in  (82)  holds  only  If  H(z)  is  stable;  that  is,  all  the  zeros  of  the  denominator 
of  (80)  must  lie  inside  O. 

In  order  to  evaluate  the  filter  ooeffiolents  j /3n } | , we  notloe  that 

wk  y*k_n  » 0 for  n > 0 (83) 

since  noise  samples  {wk}  are  unoorrelated  (see  (78))  and  filter  H(z)  is  realiz- 
able (see  (79) ).  The  first  step  we  take  is  to  express  (79)  as 


- J-  L - E fa  Vnl  • 

O L n»l  J 


Then  using  (78)  and  (83), 


wkyI“^* 

o 


Now  multiply  both  sides  of  (79)  by  y£_f  and  average;  there  follows 


upon  use  of  (81),  (83),  and  (85).  Now  if  we  let  (3n  , (86)  becomes  iden- 

tical to  (46).  Therefore,  we  can  use  solution  (63)  to  obtain  for  the  filter  ooeffi- 
olents 


= -7==rexp(-i  9),  0 £ n < p , 
>/ oo 

where  6 is  an  arbitrary  real  constant, 


TR  5303 


I 
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Substitution  oi  (87)  in  (82)  yields  the  autoregressive  speotrum  approximation 
tp  the  input  speotrum  as 


AC 


G (f)  - G (f) 
x' ' y' ' 


oo 


Ec  exp(-i2nfnA) 
no 


n=  0 


(88) 


This  is  identical  to  the  maximum  entropy  Bpeotrum  (55)  and  the  linear  predic- 
tive speotrum  (71).  The  dlsoussion  surrounding  (76)  is  relevant  here  also. 


Substitution  of  (87)  in  digital  filter  (80)  yields 


H(z) 


exp(i#) 


P 

E 


n«  0 


(89) 


This  is  stable  if  the  denominator  contains  all  its  zeros  within  O;  that  is,  H(z)  is 
stable  if  and  only  if  B(z)  of  (66)  has  no  zeros  inside  O.  This  property  has 
already  been  shown  true  in  the  discussions  under  (38),  (56),  and  (A-9). 

The  relationship  in  (86)  can  be  extended  to  j2  “ p + 1 with  the  result  that 


P 

EP  Rx1  =0,  (90) 

'll  p+l-n  ' 

n=  0 

where  Rp+i  is  now  interpreted  as  the  value  of  y^  and  was  never  speci- 

fied. If  we  combine  (90)  with  the  last  p equations  of  (86),  we  obtain 

P 

£ ^nRf-n  " °’  1 - l-  p + 1*  (91) 

n=>  0 

In  order  for  this  set  of  p + 1 linear  equations  to  possess  a nonzero  solution  for 
| 0n}p  (as  it  did  above),  we  must  have 


24 


, .•VO*Jv  .i  • 


TR  5303 


det 


R R , * • • R, 

o -1  1-p 


R,  R 
1 o 


R , R R 

p+1  p p-1 


R, 


2-p 


,0. 


(92) 


This  can  be  solved*  for  Rp+i.  But  since  this  is  identical  with  reference  19, 
equation  (1) , we  see  that  the  all-pole  digital  filter  model  is  identical  to  ohoosing 
Rp+1  suoh  that 


det 


R. 


R 


R 


p+1 


V-1 


R 


P-1 


R 


R, 


R 


-p-1 


R 


R 


(93) 


is  maximized.  Additional  interpretations  of  (93)  in  terms  of  maximum  uncer- 
tainty and  entropy  are  presented  in  references  20  and  14. 


t - ii  * 


*Of  course,  a far  more  practical  method  is  given  by  (90)  and  (87),  and 
more  generally  by  (75). 
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4.  CORRELATION  UNKNOWN;  FINITE  DATA  SET 

In  this  Beotion,  the  correlation  values  {R^}  are  unknown,  and  the  only 
information  available  about  the  random  process  x(t)  is  a finite  set  of  N samples 
Xj,  . . , xN  , from  which  we  remove  the  sample  mean.  From  these  N samples, 

; we  desire  an  estimate  of  the  spectrum  Gx(f).  Yet  we  oan  not  minimize  or  utilize 
any  ensemble  averages  as  was  done  tn  sections  2 and  3,  since  we  have  only  a 
finite  segment  of  one  member  funotlon  to  work  with. 

The  MESA  and  autoregressive  methods  of  subBeotlons  3.1  and  3.3  are  not 
easily  dlreotly  extended  to  the  case  of  unknown  correlation!  beoause  they  make 
explicit  use  of  this  correlation  Information;  see  (31)  and  (61),  respectively.  A 
direct  extension  of  these  two  methods  would  require  us  to  decide  on  the  form  of 
the  correlation  estimates  a priori,  and  oould  unnecessarily  restrict  the  quality 
of  the  speotral  estimate  we  finally  obtain.  The  linear  predlotlve  method  of  sec- 
tions 2 and  3. 2,  on  the  other  hand,  requires  that  the  ensemble  average  magnitude- 
squared  error  be  replaced  by  some  estimating  quantity  that  oan  be  readily  calculated 
from  the  available  data  {x^^  j as  a by-product,  we  may  get  estimates  of  the  cor- 
relation. Several  candidate  processing  techniques  based  on  subsection  3.  2 will 
therefore  be  considered,  and  their  processing  equations  derived.  Also,  some 
of  the  results  of  subsection  3. 1 on  MESA  will  be  adapted  and  combined  with  the 
linear  predictive  approaoh  to  form  a viable  approach  to  speotral  estimation; 
this  technique  was  originally  presented  by  Burg  in  reference  21.  In  section  6, 
all  the  techniques  will  be  oompared  by  means  of  simulation. 

4.1  YULE-WALKER  EQUATIONS 

We  begin  by  defining  In  this  subsection 

x^  ® 0 for  k<l,  k >N  , (94) 

since  these  sample,"  are  unavailable.  Taking  (58)  In  subseotion  3. 2 as  a guide,  we 
attempt  a linear  prediction  according  to 


x.  » T*  a x,  , all  k,  (95 

k « n k-n 
n=  1 

where  the  choice  of  p is  arbitrary  for  the  moment.  It  should  bo  noticed  that 
although  x^  is  defined  for  all  k,  it  is  not  expected  to  have  a good  chance  of  ac- 
curately predicting  for  k<pork>N  + 2 since  some  zero  values  of  xg  have 
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been  utilized  in  those  regions,  aooording  to  (94).  Nevertheless,  we  define  an 
instantaneous  error 


V 


p 

- £ 

n“  0 


a x.  , 
n k-n 


all  k 


<V 


1>! 


(96) 


it  is  expeoted  to  be  valid  or  meaningful,  however,  only  if  k >p  + 1 and  k <N 
(error  must  utilize  a zero  value  for  xn+i)>  Digital  filtering  operations 
(95)  and  (96)  are  stable  for  any  choice  of  coefficients  ja^  . 

Since  we  cannot  oompute  an  ensemble  average  magnitude- squared  error 
now,  on  average  magnitude-squared  error  is  defined  for  the  available  data  of 
the  single  member  function  as 


P i 

V a a*  rr  Vx,  x*  , 
m n N k-m  K-n 
m,n*0  k 


(97) 


where  £ denotes  summation  over  all  nonzero  values  of  the  Bummand  |«  ^j2  , 
regardless  of  how  meaningful  they  are.  The  normalizing  factor  1/N  is  some- 
what arbitrary;  there  are  N+p  nonzero  terms  in  the  first  sum  in  (97),  but  only 
N-p  meaningful  terms. 


We  define,  for  all  n,  m 


s = i y x,  x*  - s*  , 

n-m  N k-m  k-n  m-n 

k 


(98) 


in  which  case  (97)  yields 


F =»  T a a*  S 
„ m n n- 
m,n»  0 


m 


(99) 


This  relation  uses  Sp  only  forlil  < p.  In  order  to  minimize  F by  choice  of  filter 
coefficients  j } P , we  compute 


dF 


P 


, TT  = y S,  a , 1 < i < p 

da;  i-m  m - - ‘ 

9 m-0 


(100) 
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The  optimum  ooeffioienta  jt^|P  are  therefore  solutions  of  the  p linear  equations 
P 

£ S,  E“0,l<i<p  (a  »a  =-l)  (101) 

l-m  m ' o o ' ' 

m®  0 


or 


\ 

t 


i 

i 

l 

i 


i 


! 


P 


E 

m®  1 


S.  Z ® 
I- m m 


1 <i  < p . 


(102) 


These  are  the  Yule-  Walker  equations  for  the  optimum  filter  ooeffioienta . The 
method  here  is  called  the  autocorrelation  method  in  reference  1.  (Ab  an  aside,  in 
analogy  to  subsections  2,2  and  3.2,  identically  the  same  equations  (102)  result 
when  % is  predicted  on  the  basis  of  p future  values,  rather  than  p past  values 
as  was  done  here  in  (95);  see  (5)  and  (17)  et  seq,  and  (77)  et  seq.) 

The  minimum  value  of  average  error  F is  obtained  by  substituting  (101)  in 
(97)  and  (99); 


1 2 P P P 

| £ irj  - r k*  r s k =ir  £ s ir 

N I k n n-m  m o -m  m 

k n=  0 m=0  m=0 


P P 

£ s*  z = s - £ s*  ■ac 

m.  m o . 


m m 


(103) 


m®  0 m= 1 

where  we  have  employed  (98)  and  (101). 

Thepxpmatrix  ? on  the  left  side  of  (102)  has  the  form  of  a legal 

correlation  matrix  in  that  it  is  Hermition,  Toeplitz,  and  nonnegative  definite. 
The  last  property  follows  from 


P 

£ s, 

f,mal 


,a„a. 


= ” E 

N k 


<c~*  * 

, **  , am«f 

N £ 

f,m  = l 

k 

P 

2 

E « x, 

, m k-m 

> o 

n=  1 

(104) 
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i 


for  any  {«mf  l • Since  (104)  is  greater  than  zero  with  probability  one,  (102) 
will  possess  a solution  with  probability  one. 

A convenient  method  of  obtaining  this  solution  1b  to  combine  (101)  and  (103) 
to  get 


- £ 


m=»0 


ft 

-m  m 


- F,  l(o,  04  J <p 


Written  out  in  detail,  this  is 


s s ...  s 

“ 1 “ 

“F  1 

o -1  -p 

0 

s,  s 

-ft, 

0 

1 0 

1 

a 

• • 

• • 

• • 

• 

• 

• 

ft 

« 

ft 

s s 

-a 

0 

p 0 _ 

L pj 

_ _ 

(105) 


(106) 


(The  (p+l)  x (p+1)  matrix  in  (106)  is  nonnegative  definite,  as  a simple  extension 
of  (104)  shows.)  But  (106)  is  identical  in  form  to  (A-3),  and  the  recursive  solu- 
tion presented  in  (A-4)  through  (A-7)  applies  directly. 

The  spectral  estimate  we  adopt  follows  from  results  (68)  through  (71)  in  sub- 
section 3.2  on  linear  predictive  filtering  for  known  correlation  values:  first, 
the  optimum  transfer  function  leading  from  jxk  J to  minimum-error  sequence 
{Tk|  in  (96)  is  4 


P 

A(f)  “ £ a exp(- i2nfnA) . (107) 

n=  0 

However,  we  have  a problem  in  trying  to  accurately  estimate  the  average  mini- 
mum-error power  that  would  be  used  in  the  numerator  of  the  assumed  white 
spectrum  for  the  error  in  (70).  Although  minimum  average  error  F0  of  (103) 
could  be  used,  it  is  not  recommended  because  not  all  the  error  terms  in  the 
sum  in  definition  (97)  aro  meaningful.  Therefore,  because  of  our  inability  to 
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accurately  estimate  the  average  minimum-error  power  in  this  case,  we  shall 
adopt  as  our  spectral  estimate 


G (f)  = 
x' ' 


A'l 


A(£) 


2 £ exp(-i2irfnA) 
n=  0 


If  I < 


2A  ' 


(108) 


This  is  tantamount  to  assuming  the  average  minimum-error  power  equal  to 
unity  (in  addition  to  assuming  the  minimum-error  spectrum  white) . This  pro- 
cedure also  eliminates  level  perturbations  in  the  spectral  estimate  (108)  due  to 
random  fluctuations  in  the  absolute  level  of  the  sample  Bet  |xn|  ^ ; that  is,  from 
(102)  and  (98),  it  is  seen  that  the  optimum  values  of  the  filter  coefficients, 

{ anf  ^ , would  be  the  same  if  {Kxnj  ^ were  the  available  samples,  for  any  K. 
Therefore,  estimate  Gx(f)  in  (108)  is  also  independent  of  the  absolute  level  of 
the  available  samples.  The  choloe  (108)  allows  for  convenient  comparisons  of 
the  spectral  estimates  obtained  by  the  various  methods  presented  here. 


As  an  alternative,  (108)  could  be  numerically  integrated  over  (-^  , ^ ) , 

and  then  (108)  could  be  scaled  so  that  the  area  under  the  estimated  spectrum  is 

equal  to  the  sample  power,  4 § ixni2,  if  desired. 

^ n=l  “ 


The  implications  of  the  assumption  in  (108)  that  the  minimum-error  se- 
quence has  a white  spectrum  are  investigated  in  appendix  E . It  is  shown  that 
the  sample  crossoorrelation  between  input  sequence  and  minimum -error 

sequence  , defined  for  the  available  data  of  the  single  member  function  as 


D,  = 


- y 

N ^ 
k 


all  A 


(109) 


is  assumed  to  satisfy 


D,  =0,  1<  SL  ; 


(110) 


that  is , the  minimum-error  sequence  is  uncorrelated  (on  a single  member  func- 
tion basis)  with  all  the  past  input.  It  is  also  shown  that  the  quantities  Sf  de- 
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fined  in  (98)  (of  which  only  Sj  for  |i|  < p were  uBed  in  (99)  et  seq. ) can  be 
estimated  for  i > p + 1 according  to 


a_  S , i > P + 1 . 
n f -n 


This  relation,  (111),  which  may  not  be  true  for  the  quantities  S,  actually  ob- 
tained  from  data  |xn|  1 via  (98),  is  due  directly  to  the  assumption  that  the 
sample  spectrum  of  the  minimum-error  is  white;  see  appendix  E,  The  recur- 
sion relation  (111)  is  stable,  according  to  appendix  D,  if 


1 - £ \ z (Hi?) 

n=  1 

possesses  all  its  zeros  within  O,  But  since  matrix  in  (102)  has  the 

form  of  a legal  correlation  matrix,  we  appeal  directly  to  reference  1,  page 
567,  to  state  that  this  property  does  indeed  hold.  Therefore,  recursion  (111) 
is  stable. 

It  is  worthwhile  noting  that  no  direot  estimation  of  unknown  correlation  val- 
ues {Rk}  was  attempted  in  this  approach;  rather,  we  minimized  the  average 
error  defined  in  (97)  and  solved  directly  for  the  filter  coefficients  in  (102). 
However,  if  we  rewrite  (105)  in  the  form 


P 

£ S.  'S  - - F 
l -m  m o 

m = 0 


i < P, 


and  compare  with  (C-3),  we  see  that  the  quantity  Sf  could  be  adopted  as  an  esti- 
mate of  R|for  U2 1 < p;  that  is,  using  (98),  we  could  say 

Rf  = s<  = | £ xkx*_f  , I i I < p , (114) 

k 

(and  then  (111),  with  R replacing  S,  could  be  used  to  estimate  R#  for  I j8  I > p + 1, 
rather  than  (98)).  This  is  in  fact  the  approach  adopted  by  some  authors;  see, 
for  example,  reference  2,  equation  (19).  However,  (114)  yields  biased  esti- 
mates because 


j-  4 «•  4*- ■ t- 
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i 

| 


I 


i 


I 

> 

\ 


N-lil 


Rj  , lil  < N 
0,  otherwise 


(115) 


It  is  interesting  to  note  that  if  (114)  were  adopted  a priori  as  estimates  of 
the  unknown  correlation  values  {R^f  , then  the  MESA  and  autoregressive  ap- 
proaches of  subsections  3.1  and  3.3  oould  be  utilized  direotly,  if  the  right  sides 
of  (31)  and  (81)  were  replaced  by  . The  spectral  estimates  would  then  be 
given  by  results  identical  to  (108) , except  for  a scale  faotor.  The  major  draw- 
back of  this  approaoh  is  the  need  to  commit  oneself  to  a specific  form  for  the 
correlation  estimates , such  as  (114) , rather  than  letting  the  teohnique  itself 
yield  alternative  estimates . The  specifio  form  used  for  the  correlation  esti- 
mates oould  limit  the  quality  of  the  spectral  estimate  obtained;  this  contention 
is  proven  true  by  simulation  in  section  6. 


4.2  UNBIASED  VERSION  OF  YULE- 
WALKER  EQUATIONS 

One  method  of  obtaining  unbiased  estimates  of  the  correlation  values  {Rj| 
is  to  define  estimators 


N 

5 nTj  E Vk-I  - 5hT  „ ? Vk-( for  0 s J 5 » • <U6> 

k k«=l+l 

A rs  . 

Of  course  R_|  s R*  . These  oould  then  be  used  in  (102)  in  the  form 

P 

E **e  15J5''  <U7> 

m = l 

to  solve  for  the  filter  coefficients  { atnf  ^ . And  (108)  could  again  be  adopted  for 
the  spectral  estimate.  The  solution  for  the  coefficients  in  (117)  minimizes  no 
error  criterion;  it  merely  utilizes  unbiased  correlation  estimates.  The  dis- 
cussion under  (115)  is  relevant  to  this  approach;  how  good  the  teohnique  is  will 
be  ascertained  in  section  6. 

The  matrix  I Rj  _m]  £ of  estimated  correlation  values  on  the  loft  side  of 
(117)  is  Hermitian  and  Toeplitz;  however,  it  is  not  necessarily  nonnegative 
definite.  (This  last  property  is  shown  by  considering  the  example  p = 2,  N « 3, 
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with  xi  « 2,  x2  ■ X3  » 3;  for  then  R0  ■ 22/3,  and  R^  - 15/2. ) The  recursive 
solution  of  appendix  A could  again  be  applied  to  a modified  form  of  (117) ; see 
(106)  and  (106).  If  the  reoursive  technique  in  (111)  were  utilized  to  extrapolate 
R;  according  to 

P 

Hj  = £ i *P  + 1 <ll8> 

n®  1 

and  (116),  it  need  not  be  stable  unles  p is  nonnegative  definite.  Even  if 

(116)  were  unstable,  (108)  oould  still  be  used  as  a spectral  estimate  of  Gx(f); 
there  would,  however,  be  a greater  tendency  of  Borne  pole-pairs  of  (108)  to 
drift  olose  to  the  unit  air  ole,  O,  in  the  z-plane  and  give  rise  to  spurious  large 
peaks  in  the  speotral  estimate.  This  tendenoy  is  reduced  for  stable  recursions 
(118),  that  is,  if  (112)  possesses  all  its  zeros  within  O. 

4.3  LEAST-SQUARES  ESTIMATES  OF 
BOX  AND  JENKINS 


In  reference  16,  appendix  A7.5,  a likelihood  function  approach  to  estima- 
tion of  the  coefficients  in  an  all-pole  (that  is,  autoregressive)  filter  model  for 
generation  of  the  prooess  |xn}  is  considered.  The  end  result  (in  our  notation) 
is  given  in  (A7. 5.7)  for  real  data  by 


1 1 

Sij  2 N Di+l,J+l  = N £ VkVk’  0 “ i(  J - p 

k=  1 

and  in  (A7.5.15)  by 


(119) 


£si,VVl^p 


(120) 


This  constitutes  p linear  equations  in  the  p unknowns  |aj}  p . The  matrix  [ Sij]  P 
occurring  in  (120)  is  symmetric,  not  necessarily  Toeplitz,  and  not  necessarily 
nonnegative  definite . (The  last  property  is  shown  by  considering  the  example 
N = 5,  p » 2,  with  x2  ® x3  =>  X4  a 1,  for  then  = 3/6,  S12  = S2i  = 2/5,  S22  = 
1/5,  and  the  determinant  is  - 1/25.)  The  quantities  jSij}  also  yield  biased 
estimates  of  )R^_j[  , because 
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r ■ -N--  a,- 

1J  N 


Nevertheless  we  will  adopt  <108)  for  our  speotral  estimate  in  this  case.  The 
faot  that  we  enoounter  a non-Toeplitz  matrix  in  (120)  disallows  the  use  of  the 
recursive  technique  for  solution  in  appendix  A. 


If  the  solution  to  (120)  is  substituted  in  (112),  the  zeros  need  not  all  lie  in- 
side O,  Therefore,  there  would  be  a greater  tendenoy  for  some  pole-pairs  of 
(108)  to  drift  dose  to  O than  when  all  the  zeros  must  lie  inside  O,  as  for  sub- 
section 4.1. 

4.4  APPROXIMATE  MAXIMUM  LIKELIHOOD 
ESTIMATES  OF  BOX  AND  JENKINS 


This  teohnique  is  a slight  modification  of  the  previous  one  in  subsection 
4.3.  Namely,  in  reference  16,  under  (A7.5.18),  the  coefficients  are  solutions 
of 


- slo,  l<t<p, 


where  (see  (119) ) 


S a — D = i 

ij  ~ N-i-j  it-l,j+l  N-i-j 


N-i-j 

E 


VkVk  • 0 4 *■ 1 s p • 


The  matrix  iSij]  V ocourring  in  (122)  is  symmetric,  not  necessarily  Toeplitz, 
and  not  necessarily  nonnegative  definite.  (The  last  property  is  shown  by  con- 
sidering the  example  N *»  5,  p = 2,  with  X£  - 2,  X3  = 1,  X4  = 2,  for  then  « 3, 
S12  = S21  " 2,  S22  =>  1,  and  the  determinant  of  these  ooofficients  is  -1).  The 
quantities  {s^j}  yield  unbiased  estimates  of  {Rj_jf  {however,  every  element  in 
a particular  diagonal  can  be  different,  even  though  they  are  estimating  the  same 
quantity.  Also,  the  number  of  terms  (in  the  sum  in  (123) ) along  a particular 
diagonal  varies  with  the  position  of  the  element,  thereby  yielding  differing  de- 
grees of  stability.  Equation  (108)  can  be  used  with  (122)  to  obtain  the  speotral 
estimate.  Recursive  solution  of  (122)  is  not  allowed  because  of  the  non-Toeplitz 
character  of  the  matrix  [Sij]  £ . The  comments  at  the  end  of  subsection  4,3 
are  relevant  here  also. 
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4.D  PREDICTION  USING  VALID  ERROR 
POINTS 


The  method  of  subsection  4. 1 utilized  an  average  error  measure  defined  over 
all  nonzero  error  terms;  see  (97).  However,  as  noted  under  (96),  instantaneous 
error  < ^ is  meaningful  only  if  k > p +1  and  k £N.  Here  we  define  an  average 
magnitude-squared  error  by  summing  only  over  the  set  of  valid  error  points: 


F 


1 

N-p 


N 


E 

k-p+1 


(124) 


There  are  N-p  terms  in  this  sum.  This  procedure  is  tantamount  to  not  running 
off  the  edges  of  the  available  data  jxn}^  . Employing  (96),  (124)  oan  be  writ- 
ten as 


F 


P 


E 

m,n»0 


a a*S 
m n 


nm  ’ 


(126) 


where 


S 

nm 


1 

N-p 


„ 2 Vmxk-n 

k=p  + l 


8*  . 
mn 


(126) 


This  quantity  always  contains  N-p  terms  for  0 < n,  m < p . In  order  to 
minimize  F,  we  compute 


P 

y S,  a , 
ts  „ t m m 
m“  0 


1 < i < p . 


The  optimum  predictive  coefficients  are  therefore  solutions  of 


P 

E 


m=  0 


Se  a 
fm  m 


= 0, 


1 < i < p 


(a  = a w- 1) , 
' o o 


(127) 


(128) 


A 


I 

I 
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£ S.  a - s,  , 1 < J}  < p. 
#m  m to 


The  method  here  ia  called  the  covariance  method  in  reference  L. 

The  minimum  value  of  the  average  error  F ia  obtained  by  substituting  (128) 
in  (124)  and  (126): 

N p P P 

F a “ — £ ft  I2  - £ a*  £ S 'S-r£s  * 

o N - p . I k I n "A  ran  m o om  m 


k=p  + l 


no  0 m=0 


£ S’*  a - S - y R*  Z , 
mo  m oo  *-»  mo  m 


where  we  have  used  (126)  and  (128). 

The  p x p matrix  C S cml  i on  the  left  aide  of  (129)  is  a legal  correlation 
matrix  in  that  it  is  Hermitian  and  nonnegative  definite.  The  last  property  fol- 
lows from 

P p N 

£ S a a*  a £ a a*  rr- — £ x.  x* 

fm  m t m t N -p  , k-m  k -t 

t, m = l f,m=l  k=p  t-1 

= 2 E-.v.***  <131 

k=p  + l m«l 

for  any  Since  lS^m]  is  not  necessarily  Toeplitz,  however,  the  re- 

cursive solution  in  appendix  A is  not  applicable.  Numerical  computation  of 
[S*m]  is  eased  by  taking  advantage  of  a recursive  relation  between  Sf+lj  m +1 
and  S|m. 

The  spectral  estimate  we  adopt  is  given  by  (108),  However,  note  that  we 
could,  if  desired,  get  an  estimate  here  of  the  average  minimum  error  power 
E0,  used  in  (70),  according  to 


...  ^ \ . ' 1 ’ flA"***  ''Z1  ’'C*  * 
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This  quantity  is  meaningful  because  (130)  involves  only  the  valid  error  terms. 

Equation  (129)  is  similar  to,  but  not  ldentioal  with,  the  form  of  (117).  The 
quantities  {8^}  . defined  in  (126),  yield  unbiased  estimates  of  {Rn-mf  ; how- 
ever, 9very  element  in  a particular  diagonal  can  be  different,  even  though  they 
are  estimating  the  same  quantity. 

If  the  solution  to  (129)  is  substituted  in  (112),  the  zeroB  need  not  lie  inside 
O,  despite  the  nonnegative  definite  property  demonstrated  in  (131).  (The  ex- 
ample 


p = 1,  N » 2,  yields  a^  «*  X2^xi  (133) 

and  gives  a zero  location  of  (112)  whloh  can  lie  anywhere  in  the  z-plane.) 
Therefore,  the  comments  at  the  end  of  subsection  4.3  are  relevant  here  bIbo. 

4.6  FORWARD  AND  BACKWARD  PRE- 
DICTION USING  VALID  ERROR  POINTS 

It  was  noted  in  subsections  2.2,  3.2,  and  4.1  that  "prediction"  based  on 
future  values  of  the  input  {x^  | yielded  an  equivalent  spectral  estimate  to  that 
obtained  by  prediction  based  on  past  values.  Here  we  oombine  the  two  tech- 
niques . The  forward-predicted  value  of  xk  is 

P 

\ s £ afi  xk_n , p + 1 £ k £ N , (134) 

n=  1 

where  we  limit  k to  the  range  [p  -i- 1,  N) , in  anticipation  of  the  fact  that  we  can 
only  measure  valid  errors  in  that  range;  see  (96)  et  seq.  The  backward-pre- 
dicted value  of  x^  is 


P 

*kB  E ■;  Vn’  ^k*N-P'  <135> 

n=  1 

where  we  again  limit  the  range  of  k.  (See,  for  example,  (15),  (22),  and  (77).) 
The  forward  and  backward  errors  are,  respectively, 
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W \ ‘ £ an  Vn  ’ P + 1^k^N- 
n=  0 

(136) 

Tk  " *k  - Xk  “ £ ‘JVn-  l^k^N-p. 
n«  0 


where  a ■ - 1. 

o 


An  overall  average  magnitude- Bquared  error  Is  defined  aa 

N-p 


F s 


N 

E 

k«p  + l 


w=$ 

where,  In  this  subsection, 
1 


k 2 

k 


- 2 

k 


>r)=  £, 


a a*  8 , (137) 

m n nm  ' 


N N-p 


S 


nm  2 (N-p) 


„ 2 , Wk-n  +£  Wl 

k=p+l  k-1 


* 

k+ml  * 


(138) 


This  quantity  always  contains  2 (N-p)  terms  for  0 < n,  m < p . Two  useful  prop- 

(139) 


ertles  of  are  immediately  available: 


s » s*  , s - s* 

mn  nm  p-n,  p-m  nm 


These  properties,  plus  a recursive  relation  relating  Sn+if  m+i  and  Snm , ease 
the  numerical  computation  of  matrix  [ 8nm]  . 

We  minimize  F by  choioe  of  jam^  , getting  (see  (127)  - (129) ) 

P 

E & * a Kn  * 1 * n * P • <140> 

M nm  m no 
m = 1 


The  minimum  value  of  F is  (see  (130) ) 
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F « S - S*  tf  . (1- 

o oo  mo  m ' 

m«l 

The  method  here  Is  an  extended  version  of  the  oovarianoe  approach  in  refer- 
ence 1. 

The  matrix  [Sfm]  ^ is  Hermltlan  and  nonnegative  definite: 


P . / N p 2 N-p  P 2\ 

2 s<mV*l  ":2(N-p)(  £ £ “mXk-m  + £ ^ am  k+m  )i( 

f,m-l  \k=p  + l m=»l  k-1  m-1  ^14J 

for  any  . However,  this  matrix  Is  not  necessarily  Toeplltz;  therefore, 

we  oannot  apply  the  recursive  solution  of  appendix  A. 

i1  1 * 

The  speotral  estimate  we  adopt  is  obtained  by  substituting  the  solution  of  . 
(140)  in  (10B).  An  estimate  of  the  average  minimum  error  power  Eol  used  in 
(70),  is  available  here  according  to 

E » F , (14! 

oo  ' 

If  desired,  where  FQ  is  given  by  (141).  This  is  meaningful  because  (137)  util- 
ized only  the  valid  error  terms. 

Ir.  analogy  to  (126),  the  quantities  jsnm|  In  (138)  yield  unbiased  estimates 
of  {Rn-m}  • Nevertheless,  if  the  solution  to  (140)  is  substituted  in  (112),  the 
zeros  need  not  lie  within  O,  despite  the  nonnegative  definite  property  shown  in 
(142).  For  p ■ 1,  we  find 


Vl  +X3X2  + •••  + Vn-1 
” L 12  . I 12  .1 


11  2 ^ 2 2 

2 X1  + X2  + X3  + •**  + XN-2 


I 12  11  12  * 


And  sinoe 
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It  follows  that 


|^|  < 1.  (HQ) 

So  for  p ■ 1,  the  zero  of  (112)  must  He  within  O,  (unless  x^  ■ A exp(ikB)  for  all 
k,  in  whioh  case  it  lies  on  O).  However  for  p « 2,  N ■ 3,  and  real  data,  the 

zeros  of  (112)  lie  at  r + Jr  - 1,  where  r ■ — g . So  if  Irl  £ 1,  both  zeros 

lies  on  O,  whereas  if  Irl  >1,  one  zero  lies  outside  O.  Therefore,  the  oom- 
ments  at  the  end  of  subseotion  4.3  are  relevant  here  also. 


4.7  BURG  TECHNIQUE 


The  key  to  this  technique,  first  presented  in  reference  21,  is  the  observa- 
tion from  equation  (A-6)  in  appendix  A that  if  the  particular  p-th  order  coeffi- 
cient afr)  op  be  evaluated,  the  rest  of  the  p-th  order  predictive  filter  coeffi- 
cients, ajjjr  , 1 <,  k£  p - t,  could  be  evaluated  from  (p  - l)-th  order 
coefficients , This  relation  (A-6)  holds  true  for  the  solution  of  normal  equations 
^-3)  even  if  iR^}  are  replaced  by  estimated  values.  Explicitly,  if  estimates 
R0,  ftj_,  ...,  Rp-i,  and  ajP)  are  considered  known  in  the  matrix  equation 


R £ 1 £ 

i 

"p(pT 

o -1  -p 

A A 

R.  R 

-a» 

0 

1 o 

• t 

1 

• 

zz 

• 

• * 
a • 

& & 

• 

• 

-a<» 

4 

0 

L p oJ 

L pj 

H>  mm 

(147) 


then  we  have  p + 1 linear  equations 
P^.  (Notice  that  whereas  Rp  was 
tion  is  reversed  here  for  these  two 
P >1*  by 


in  the  p + 1 unknowns  a, 
<P)  1 


(P) 


JP)  £ 

**p— 1*  "p* 


known  and  unknown  in  (A-3),  the  situa- 
variables.)  The  solutions  are  given,  for 


,(P-D< 

p-k 


k 


1,  2,  ....  p - 1 (no  terms  if  p=  1)  (148) 
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R - £ R . a^ , 
p p-k  k * 

p - 

P 

L 

S a^ 

Vk  k * 

(149) 

k-1 

k-1 

p(P)  . p(P-D 

(- 

P 

y 

(150) 

The  quantities  |pk|  in  (149)  are  the  estimated  normalized  correlation  ooeffi-  r>( 
oients  |Rk/R0|  . The  recursion  (150)  is  started  with 


N 2 

K - N £ |Xnl  • 


n-1 


(l&l) 


which  is  the  sample  power  of  the  available  samples.  ' A method  of  evaluating 
for  p £ 1 is  treated  below. 

The  method  presented  here  is  a combination  of  references  21  and  7.  It 
begins  by  defining  zero-th  order  forward  and  backward  sequences  according  to 

fi0)  “ xn  • bf-x  . l^n  . (152) 

n n n n 

The  p-th  order  forward  and  backward  sequences  (residuals)  for  p ;>  1 are  de- 
fined according  to 


»« a fi-1'  - g 

n n °p  n-1 

bw  - b(p:1)  - s* 
n n-i  °p  n 


for  p + 1 <;  n < N . 


(153) 


(These  can  be  Interpreted  as  one-step  forward  and  backward  prediction  errors.) 
A chain  interpretation  of  (153)  is  presented  in  figure  3.  (From  the  known 
correlation  results  in  subsection  3,  2,  if  we  define 


x - V)  ar 

n lfcl  * 


(P) 


n-k 


we  find  that  figure  3 results,  with  g 


replaoed  by  a^. ) 


...» 
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Figure  3.  Chain  Interpretation  of  Burg  Technique 

The  average  magnltude-aquared  value  of  the  p-th  order  forward  and  back- 
ward sequences  Is 


u-  p + 1 


L_  £ /Lip-d  „ t(p -i)  2 ^ h(p-u  „*  f<p-u|2'\ 

2(5h3  E \|f„  -S„bn-1  * \-l  'gpfn  >pal- 

n«p  + l (15 

We  wish  to  minimize  this  average  power  at  the  p-th  stage  by  choice  of  cross- 
gain  gp,  We  find  the  optimum  choice  to  be 


2 £ 

n°p  + l n 11 

E i (j^W) 

n«p+l  ' 


Num(p) 

Den(p) 
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When  (155)  is  substituted  in  (153),  the  results  are  called  the  residuals.  The 
minimum  value  of  the  residual  power  at  the  p-th  stage  is  obtained  by  substitu- 
ting (155)  in  (154)  and  is  expressible  as 


i: 


M 


(156) 


The  quantities  neoessary  for  this  evaluation  are  available  when  (155)  is  evalu- 
ated. The  value  of  (156)  will  never  be  smaller  than  (141),  slnoe  the  procedure 
here  1b  a step-by-step  procedure,  not  a simultaneous  procedure  as  used  in 
subsection  4. 6. 


An  Immediate  recursion  tor  the  transfer  funotions  of  the  p-th  stage  in  fig- 
ure 3 is 


3 (l>)(z) , 2F'p-i;(Z)  - 8l-ir'(!) 
B^z)  - a“le*"l><z>  - g*g(p“l)(z) 


f<p-l), 


with  starting  values 


If  we  let  transfer  functions 


3(0)(/,)  - ®(0)(n)  - 1 . 


(P)  “k 


3 W«  - 1 - £ a-  ' z 
k33 1 


» P S;  1 . 


(157) 


(168) 


P-1 


B^W  - - E «®k*  " s'^VV 


k»  0 


(159) 


the  solution  is 


„(P)  M „ V . 

ttp  " v p 


(160) 


with  the  lower  order  coefficients  given  by  (148).  Thus,  the  only  remaining 
quantity,  ajf' , that  web  necessary  for  solution  of  (147)  - (150)  is  given  by  (160) 
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and  <155),  along  with  (152)  and  (153).  To  the  three  lowest  orders,  the  solutions 
are  given  by 


N 


p = o,  P(0)  « ft  « £ £ ix 
r o N *-*  n 

n=l 


p = 1,  a 


(1) 


A A (1) 

R,  = R aj  ' 
1 o 1 


p(D  = p<°> 


P “ a.  ai2)  = g, 


(>  - KT) 


2 2 

(2)  ~ a(1)  - a(2)  a(1)* 
1 1 2 1 


R a(2)  + ft  a(2) 

z R1  ai  + R0  a2 

pW.pW  (l-|af|2). 


(101) 


(162) 


(163) 


It  will  be  ohserved  that  for  p = 1,  is  identical  to  (144);  In  fact,  the  proce- 
dures are  identical  in  this  case.  It  should  also  be  noted  +hat  at  each  stage,  an 
estimate,  Rp,  of  the  true  correlation  value  Rp  becomes  available  via  (149),  and 
is  unchanged  by  the  addition  of  any  further  stages  (larger  p) . 

It  was  demonstrated  in  (A-  9)  that  the  magnitude  of  was  bounded  by 
unity  if  the  known  correlation  matrix  R was  nonnegative  definite.  The  same 
property, 


(164) 


is  true  here  in  the  case  of  unknown  correlation  when  r^  is  determined  by  (160) 
and  (155);  see  appendix  F.  This  is  sufficient  to  show  that  all  zeros  of  (A- 10) 
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lie  inside  6;  see  reference  11,  for  example.  Therefore,  the  recursion  (149) 
can  be  used  in  the  form 

ft,  ■ h ft,.*  . S > P + 1 , 

k«l 

to  extrapolate  the  estimated  correlation  values  beyond  p-th  order,  with  the  p-th 
order  coefficients  {ajjP^f  ^ , and  is  guaranteed  to  be  stable.  Division  of  (165)  by 
R0  yields  the  normalized  correlation  coefficients.  Recursion  (166)  is  similar 
in  form  to  those  encountered  in  (75),  (111),  and  (118). 

The  quantity  P^  that  results  as  the  solution  (150)  of  matrix  equation  (147) 
is  not  the  minimum  average  magnitude- squared  error  as  it  was  for  known  cor- 
relations see  (A-3) , (A-7),  and  (67).  In  fact,  pO?)  has  no  direct  physical  sig- 
nificance; it  is  merely  the  variable  left  over  in  that  position  in  the  normal  equa- 
tions (147)  when  modified  from  the  case  of  known  correlation  values,  (A-3). 
Rather,  F<f  in  (154)  and  (156)  is  the  minimum  average  magnitude- squared 
error  of  the  forward  and  backward  residuals,  (153),  of  the  available  data. 

Thus,  as  far  as  picking  an  "optimum"  value  of  p at  which  to  terminate  the  re- 
cursion in  (147)  - (150)  is  coAoerned,  the  latter  quantity  has  more  physical  sig- 
nificance. However,  the  two  quantities  are  very  dose  to  each  other  for  no  bad 
data  points,  especially  for  N-p  large;  see  appendix  Q.  Both  quantities  are 
readily  oaloulated  at  any  stafce  via  (150)  and  (156),  respectively. 

The  transfer  functions  from  input  x to  the  p-th  order  residuals  are  given  in 
(159).  Therefore,  the  speotra  of  the  residuals  are  given  by 

G^f)  - G^(f)  « | 0(p)(exp(i2TrfA))|2  Gx(f) . (166) 

Now  if  the  ohain  in  figure  3 has  been  oarried  to  the  stage  where  further  values 
of  oross-galn  g_  would  be  substantially  zero,  then  the  residuals  are  approxi- 
mately white.  Therefore,  an  estimate  of  the  input  spectrum  is  available  from 
(166)  and  (159)  aooording  to 
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where  the  residual  power  has  been  set  at  unity;  see  the  disousBion  under  (108) . 
Two  alternatives  to  this  scale  faotor  are  discussed  in  appendix  H;  namely,  it  is 
shown  that  P^P)  and  F$?)  are  both  meaningful  soale  factors  that  could  be  applied 
to  (167). 

The  estimated  correlation  values  in  (149)  are  generally  biased.  This  may 
be  anticipated  from  the  complicated  forms  of  (149),  (148),  and  (186),  shme 
additional  statistics  than  simply  need  to  be  known  in  order  that  Rp  be 

capable  of  evaluation;  that  is,  Hp  depends  on  muoh  more  than  Just  x^+pXfc , for 
the  Burg  method.  This  biasedness  may  be  proven  for  a simple  example  with 
p»l,  N <=  3.  (ft0  in  (151)  is  unbiased;  and  for  p « 1,  N « 2,  we  find  - XgxJ  , 
which  is  unbiased.)  For  real  data,  with  random  variables  being  zero- 

mean  unit-variance  Gaussian,  and  koki  ™ ™ + i,  R3XT  •»  0,  we  find  (in 


append  ix 
zero. 


t-variance  Gaussian,  andlf^T  ™ *3*2  “ ±y|«»  *»  0,  we  find  (in 

I)  that  Si  “ + = -Z  — “ + - (.9484).  The  bias  is  slight  but  no 

"v2  ® 


In  summary,  the  Burg  algorithm  for  data  processing  oonslsts  of  initialisa- 
tion (152);  followed  by  the  oross-gain  calculation,  in  (155);  filter  coefficients  via 
(160)  and  (148);  and  normalized  correlation  coefficients  (146)  (If  desired) 
at  every  stage.  The  update  required  at  each  stage  la  given  by  (163),  and 
the  extrapolated  normalized  correlation  coefficients  at  any  stage  are  available 
from  (165), upon  division  by  ftQ. 

4.8  SUMMARY  OF  PROPERTIES  OF 
TECHNIQUES 

The  solution  for  the  filter  coefficients  in  the  techniques  considered  above 
can  be  put  in  the  form 


r(-»)  - Fq  » . 


The  properties  of  the  estimated  correlation  matrix  R (if  desired)  are  tabulated 
in  table  1.  (Actually,  several  of  the  "No"  entries  should  be  "Not  Necessarily,") 

It  will  be  seen  that  none  of  the  techniques  possessor  a "Yes"  for  all  the 
properties.  The  Yule-Walker  and  Burg  techniques  possess  everything  but  the 
unbiased  property;  however,  the  unbiased  property,  per  se,  of  the  correlation 
estimates  is  not  necessarily  a desirable  feature  for  spectral  estimation,  as  will 
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Table  1.  Properties  of  Estimated  Correlation  Matrices 


Technique 

Correlation 

Estimates 

Unbiased 

Hermitian 

Toeplitz 

Nonneg- 

ative 

Definite 

Stable 

Recursion 

' 

Yule-Walker 

(114) 

No 

Yes 

Yes 

Yes 

Yes 

Unbiased 

Yule-Walker 

(116) 

Yes 

Yes 

Yes 

No 

No 

Least  Squares 
of  Box  and 
Jenkins 

(119) 

No 

Ycb 

No 

No 

No 

Approximate 
maximum 
likelihood  of 
Box  and 
Jenkins 

(123) 

Yes 

Yes 

No 

No 

No 

Prediction 

(126) 

i .... 

Yes 

Yes 

No 

Yes 

No 

Forward  and 

Backward 

Prediction 

(138) 

Yes 

Yes 

No 

Yes 

No 

Burg 

(149) 

No 

Yes 

Yes 

Yes 

Yes 

be  seen  by  later  simulation  results.  On  the  other  hand,  simultaneous  satisfac- 
tion of  the  three  properties  of  Hermitian,  Toeplitz , and  nonnegative  definite 
guarantees  that  a stable  recursion  and  nonspiky  spectral  estimates  result;  see 
reference  1,  page  567. 
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6.  CORRELATION  UNKNOWN:  FINITE  DATA  SET 
WITH  BAD  DATA  POINTS 

In  some  applications,  some  data  values  can  be  bad  as  a result  of  malfunc  - 
\ tlonlng  equipment  or  human  errors  in  reading  or  recording,  for  example.  Also, 

some  data  values  oan  be  missing  as  a result  of  equipment  being  Inadvertently 
or  intermittently  turned  off  for  calibration  purposes,  for  example;  or  some  sec- 
tions of  data  oan  be  contaminated  by  strong  burst-like  noise  and  be  virtually 
useless  in  those  sections.  All  of  these  problems  oan  be  characterized  mathe- 
matically by  saying  that  of  the  available  data  set  |xn}  ^ , the  values  xn  for  the 
distinct  integers 

n - Mx,  Ma,  ...  Mfl  (169) 


are  known  to  be  bad  (or  missing).  The  B bad  locations  {Mj}  ® are  presumed 
to  be  known.  The  bad  data  points  can  be  regularly  spaoed,  or  randomly  spaced, 
or  a combination,  depending  on  the  application,  it  will  make  no  difference  to 
the  techniques  to  be  developed  here. 

In  this  section,  we  wish  to  estimate  the  input  spectrum  despite  the  presenoe 
of  known  bad  points.  The  last  two  methods  in  subsections  4. 6 and  4. 7 will  be 
extended  to  oover  this  oase.  The  reason  we  do  not  extend  the  other  methods  in 
seotion  4 will  become  dear  when  we  compare  the  various  techniques  by  simu- 
lation in  section  6, 

5.1  FORWARD  AND  BACKWARD  PRE- 
DICTION USING  VALID  ERROR  POINTS 

The  method  to  be  presented  here  1b  very  Bimilar  to  that  given  earlier  in 
subsection  4.6;  accordingly  the  treatment  will  be  briefer.  For  convenience  and 
to  enable  a better  estimation  of  the  true  speotrum  with  a limited  number,  p,  of 
parameters,  we  subtraot  the  sample  mean  of  the  N-B  good  data  points  so  that 

N 

j^B  S *n  " ° ’ <170A) 

n-1 

where  $ denotes  that  we  skip  those  values  of  n in  the  set  (169);  that  is,  we 
simply  ignore  the  bad  data  pointB  — this  is,  in  fact,  the  main  theme  of  the 
methods  to  be  presented.  We  attempt  no  interpolation  on  the  bad  points,  nor  do 
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we  set  them  equal  to  zero  or  the  sample  mean.  We  also  scale  the  good  points 
so  that  the  sample  variance  is  unity: 


1 

N-B-l 


N 


S 

n=  1 


(170B) 


This  helps  avoid  overflow  and  underflow  problems  in  the  numerloal  manipulation 
of  large  arrays  encountered  for  large  p. 

A forward  predlotion  of  is  afforded  by 

P 

E an  *k_n  , P + 1 < k < N , (171) 

n®  1 

provided  that  k-1,  k-2 k-p  4 Mj_,  M2,  ...»  Mg.  Then  a valid  forward 

error  oan  be  defined  as 


tksSk'\*  T anVn  (“o“ -l)iP  + lSk<N,  ,172) 
n“  0 

provided  that  k,  k-1,  ....  k-p  4 M^,  Mg,  ...»  Mg  ; that  is,  is  defined  for 
p + 1 < k < N exoept  for  k in  the  set  of  integers 

Ml’  M^  + 1,  . . . , M + p 

Mg , Mg  + 1 , ...  i Mg  + p 


Mb.  Mb+1,  ....  . M0  +p  . 

If  any  numbers  in  set  Ip  are  <p  + l or  > N,  they  axe  not  enoountered  in  the  error 
definition  (172).  Let  Bp  denote  the  number  of  distinct  integers  in  L whioh  are 
;>  p + 1 and  <N;  this  is  the  number  of  gaps  (bad  points)  in  the  error  soquenoe 
(172). 
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We  now  define  an  average  forward  error  over  the  valid  error  points  aa 


A 

F 


1 £ .2 

N-p-B  , ^ < |k| 

P k=p  + l 


(174) 


where  $ denotes  "not  contained  in, " and  N -p  - Bp  is  the  number  of  terms  in  the 
sum.  Substituting  (172)  in  (174),  we  obtain 


^ ^ aman  N-p-B  ^ Xk-m  Xk-n  ’ 

m,n=  0 p k=p  + l 

kf'P 


(176) 


A backward  prediction  of  is  available  as 

P 

V £ <Vn-  1Sk£N  ’p'  (170) 

n=  1 

provided  that  k+1,  k+2 k+p  ^ Mi,M2,  ...»  M0.  And  a backward  error 

W\-  £ a;v«  V-1''  <”7> 

n=>  0 

is  available  if  k,  k+1,  . . . , k+p  ^ M],,  M2>  . . . , Mg.  Letting  Jl  =■  k+p  in  (177), 
we  can  write 


« = V 

i-p  " 
n=  0 


X A , 

n i-p+n 


p + 1 < .&  < N , 


(178) 


if  i is  not  contained  in  the  set  Ip  defined  in  (173).  The  we  can  define  an  average 
backward  error  over  the  valid  error  points  as 
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v 

F 


N-p-B 


N 

E 

p i-p+i 
i 


■p+] 

fh> 


2 


P 

- E 

m,n»  0 


* 1 

a a — — — 

m n N-p-B 

* w 


N 

E 

i»p+i 


x,  , XT 
l-p+n  f • 


p+m 


(179) 


(180) 


where  we  have  substituted  (178). 


We  are  now  in  a position  to  define  an  overall  average  error  as 


F - J (F  + ¥)  - £ a a*  S 

2 m n nm 

m,n*  0 


(181) 


where,  from  (175)  and  (180), 


N 


S ra  uni  ■ .■ 

nm  2(N 


- p - B ) ^ (Xk-m  \-n  + Xk-p+n  ^-p+m) ' 

p k“p+l  ' ' 

kflp 


(182) 


It  should  be  noticed  that  (182)  does  not  tell  us  merely  to  sum  over  the  "good" 
products,  but  rather  to  exclude  set  Ip.  The  number  of  terms  in  the  sum  (182) 
is  the  Home  for  0<n,  msp  and  is  N-p-Bp,  (For  no  bad  points,  (182)  reduces  to 
(138) . ) Two  useful  properties  of  Snm  are 


S 


mn 


S*  , 8 

nm  p-n , p-m 


S*  . 
nm 


(183) 


The  quantity  Snm  is  aft  unbiased  estimate  of  Rn>.m;  however,  the  presence  of 
bad  points  will  increase  the  variance  of  Snm;  see  reference  5.  The  matrix 
t snmJ  i 1®  Kermltian  and  nonnegative  definite. 

The  optimum  predictive  filter  coefficients  {'Em}  ? are  obtained  by  mini- 
mizing (181): 
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T]  S a = S , 1 < J2  < p . 
t mm  in 


m = 1 

The  minimum  average  error  is 


to 


(184) 


P 

£ 

m=  1 


F = S - J]  S*  a 
o oo  ^ mo  m 


(186) 


And  since  the  sample  variance  of  the  good  data  points  was  set  equal  to  unity  in 
(170B),  (186)  is  a relative  error  measure  that  can  be  uBed  to  decide  what  value 
of  p should  be  used  in  (171)  and  (176);  see  reference  1,  equations  (41)  and  (89) 
et  seq.  The  spectral  estimate  we  adopt  is  obtained  by  substituting  the  solution 
of  (184)  into  (108),  as  usual.  The  quantity  F0  in  (186)  oould  be  used  as  a soale 
factor,  if  desired,  according  to 


E - F . 
o o 


(186) 


6.2  BURG  TECHNIQUE 

The  problem  setting  is  the  same  as  that  for  the  previous  subsection,  in- 
cluding (169)  - (170).  The  solution  is  identical  to  that  for  subsection  4.7,  up  to 
(160).  Now  we  define  zero-th  order  forward  and  backward  sequences  as 


f(0)  » x , b(0)  =>  x , 1 £ n £ N , n £ I 
n n n n to 


(187) 


where  we  again  employ  the  definition  (173).  The  first-order  sequences  are  de- 
fined as 


f(1)  - f(0)  - g,  b<°> 
n n el  n-1 


- if.  - B;  if 


for  2 n ^ N , n ^ I , 


(188) 


n-1  °1  n 


where  the  restriction  of  set  1^  is  due  to  the  fact  that  the  first-order  sequences 
cannot  be  formed  (evaluated)  in  set  I^.  We  choose  cross  gain  g^  to  minimize 
the  average  error 
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rw  - 1 v (|fa> 

"2<N-1-B1)  n-2^U 

n$Xl 


where  is  the  number  of  terms  in  the  sum.  The  solution  iB  given  by 


2 £ fWbw; 

nf2  n 

n^i 


1 £ (lit  • 1*1’) 


f 


With  this  value  of  g^,  we  can  now  compute  values  for  residuals  f^,  b£^  In 
(188)  and  oontinue  the  prooedure. 


At  stage  p , we  have 


V . ,r , . ^ 

b<*»  = b(p:l) . g*  t(p-l> 

n n-l  °p  n 


p + l<n<N,n^I 


The  choice  of  cross-gain  gp  that  minimizes  average  error 


f<p)  5 w-t--b  i £ rf  * i bf f)  (p(o)  - 

P n°p  + V ' 

nf Ip 
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Is 


N 

2 L 

n-p 
n 


n-p  + 1 

1 

E_ 


f(p-l)  b(p-l)* 
n n-1 


N /l  ^ — 1)|2  + |b(p-i)'2 

i+i  " 1 • n_l 


£ (|*r 

n»p+l  ' 

n|lp 


Num(p) 
Den(p)  * 


and  the  minimum  value  of  (192)  can  be  expressed  as 
■pto)  /i 


1 |8v  JBhM. 

(f(0>  - l) 

|gp|  ) 2(N-p  - Bp) 

\ o / 

(193) 


(194) 


This  is  also  a relative  erro.',  due  to  the  normalization  (170B),  and  oan  be  used 
as  an  indicator  when  to  terminate  the  recursion  procedure  in  (191). 


It  may  be  seen  from  (192)  and  (193)  that  the  sumB  are  merely  taken  over 
those  values  of  n where  the  summands  are  defined.  The  number  of  terms  in  all 
the  sums  is  N-p-Bp. 

As  in  subsection  4.7,  the  filter  coefficients  are  given  by 


a 


(P) 


V 


p > 1 


and  for  p > 2 , by 


a(P>  tt(P-l)* 
P P-k 


l <,  k < p - 1 . 


(195) 


(196) 


Equations  (147)  through  (150)  still  hold  true.  The  starting  value  of  P(0)  ia  now 
1,  by  virtue  of  normalization  (170B).  Recursion  (165)  for  J i >p  + l is  still  valid 
and  is  stable  since 


.<P> 


LSI, 


(197) 


as  may  be  seen  from  (193)  and  appendix  F.  The  speotral  estlmato  is  again 
given  by  (167).  The  discussions  in  appendixes  G and  H are  relevant  here  also. 
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6.  COMPARISONS 

All  the  techniques  considered  In  section  4 will  now  be  compared  in  terms  of 
their  resolution  capability,  bias,  and  statistical  stability,  by  means  of  a simu- 
lation approach.  In  particular,  the  fourth-order  autoregressive  process  which 
was  intensively  investigated  in  reference  2 (see  figures  4a  and  5a)  will  be  the 
basio  prooess  of  interest  here  also.  It  is  characterized  by 

4 

*k”  £ Vk-n+V  (196) 

n=  1 

0 

where 


a,  = 2.7607,  a =>  - 3.8106,  «0  = 2.6536,  a - - O.Di.38,  (199) 

1 « 3 4 

and  where  jw^}  is  Gaussian  white  noise.  We  restrict  consideration  to  real 
processes  here.  We  will  not  address  the  problem  of  how  best  to  pick  the  value 
of  p used  in  the  techniques  of  seotions4and5,  but  shall  instead  set  p equal  to 
the  known  value,  4,  and  concentrate  on  the  ability  of  the  various  techniques  to 
estimate  the  parameters  in  (199;,,  and  thereby  the  spectrum  of  jx^}  , from  a 
finite  set  of  N data  points. 

The  simulation  method  consists  of  the  generation  of  100  independent  reali- 
zations of  the  prooess  in  (198)  in  steady  state.  The  coefficients  in  (198)  are 
estimated  for  eaoh  of  the  100  realizations,  and  the  corresponding  100  estimated 
spectra  aro  computed  by  means  of  (108),  for  every  technique  in  sections  4 and 
5.  The  examples  to  be  considered  are  summarized  in  table  2,  where  N is  the 
number  of  data  points  in  each  realization  (trial),  and  B is  the  number  of  bud 
points  in  eaoh  realization.  The  corresponding  figures  are  collected  together  at 
the  end  of  this  section. 

6.1  NO  BAD  DATA  POINTS 

In  figure  4A,  the  100  different  estimated  spectra,  one  for  each  of  the  100 
independent  trials,  are  plotted  for  the  Yule- Walker  approach,  and  for  N * 40 
data  points.  In  figure  4B,  the  (power)  average  spectrum  of  the  100  estimated 
spectra  is  plotted,  along  with  the  true  spectrum  of  prooess  (198)  and  (199). 

The  true  speotrum  is  scaled  so  that.  Its  area  is  equal  to  that  of  the  average 
spectrum . It  will  be  seen  from  figure  4A  that  there  is  a great  don’  of  variabil- 
ity in  the  individual  spootral  estimates.  From  figure  4B,  we  observe  that  the 
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Table  2.  Simulation  Examples 


Figure  - Number  of  Number  ol 

Number  Data  Points  Bad  Points 

N B 


Yule-Walker 


Yule- Walker,  Unbiased 
Least  Squares  of  Box  and  Jenkins 


Approximate  Maximum  Likelihood 
of  Box  and  Jenkins 


Prediction,  Valid  Error  Points 


Forward  & Backward  Prediction 


Burg,  Uniform  NoiBe 


Forward  & Backward  Prediction 


Forward  & Backward  Prediction 


Forward  b Backward  Prediction 


Forward  & Backward  Prediction 


Forward  & Backward  Prediction 
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average  spectrum  does  not  resolve  the  two  narrowband  peaks  of  the  true  speo- 
trum*  i In  faot,  this  same  conclusion  is  true  for  the  Individual  spectra  in  figure 
4A.  A severe  bias  exists  in  the  skirts  of  the  average,  spectrum,  which  gives  a 
gross  overestimate  of  the  power  in  bardB  away  from  the  peaks.  Thus,  the 
Yule- Walker  approach  has  poor  resolution,  severe  bias,  and  substantial  vari- 
ability. 

The  corresponding  results  for  the  unbiased  version  of  the  Yule- Walker  ap- 
proach are  displayed  in  figure  5.  Rather  than  improving  the  situation,  it  is 
found  that  the  spectral  estimates  are  worse  in  every  regard.  The  spectral 
estimates  with  strong  spikes  near  f « -g-j-  are  manifestations  of  pole-pair  loca- 
tions i.'  ostimate  (108)  that  are  very  olose  to  the  unit  oirole  O.  Recall  from 
subseotion  4 .2  that  the  zeros  of  (112)  need  not  lie  inside  O;  see  the  discussion 
below  (118). 

The  unblasod  correlation  estimates  utilized  above  in  the  normal  equations 
are  of  the  same  form  as  those  suggested  in  reference  5 for  missing  data,  when 
spectral  estimation  is  attempted  directly  via  (2).  But  since  the  performance  of 
these  unbiased  correlation  estimates  is  so  poor  here,  they  aro  not  considered 
worthwhile  in  the  presence  of  bad  data  pointB,  when  spectral  estimation  is  ac- 
complished via  (108).  Whether  they  are  worthwhile  for  use  in  (2)  is  not  known. 

Results  for  the  least-squares  approaoh  of  Bex  and  Jenkins  are  given  in  fig- 
ure 6,  The  variability  is  less  than  that  for  the  Yule-Walker  estimates  in  figure 
4A.  And  some  resolution  is  aoiiieved  in  figure  6B,  in  addition  to  good  skirt 
selectivity . There  Is  still,  however,  a large  number  of  spiky  spectral  esti- 
mates, ns  anticipated  in  the  discussion  under  (121). 

Conditions  are  not  much  improved  for  the  approximate  maximum  likelihood 
method  of  Box  and  Jenkins  presented  in  figure  V.  There  happens  to  bo  one  par- 
ticular spectral  estimate  with  a very  large  spike  (a  zero  virtually  on  O)  that 
severely  influences  tho  average  power  level . The  variability  in  the  estimated 
skirt  level  is  quite  small  for  this  technique  (a<i  well  as  for  the  previous  one). 

in  figure  8,  tho  results  for  prediction  using  only  the  valid  error  points  are 
presented.  Tho  resolution  and  bias  in  figure  8B  are  observed  to  be  very  good, 
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*Thls  same  conclusion  is  reached  in  reference  2,  figure  5b,  for  the  same 
number  of  data,  points.  Increasing  p (above  4)  does  recover  some  of  the  resolu- 
tion of  the  two  narrowband  peaks,  but  it  does  not  reduce  the  severe  bias  of  tho 
Yule- Walker  approaoh. 
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and  except  for  some  spiky  estimates  in  figure  8A,  the  variability  of  the  individ- 
ual estimates  is  fairly  small. 

The  situation  is  still  better  when  we  consider  forward  and  backward  predic- 
tion, using  only  the  valid  error  points,  in  figure  9.  There  are  a couple  of  spiky 
estimates,  but  they  are  not  excessively  large,  as  they  were  previously.  The 
bias  and  resolution  are  very  good  in  figure  9B.  Although  the  zeros  of  (112)  need 
not  remain  inside  O for  this  technique,  it  was  found  that  in  all  100  trials,  no 
zeros  were  ever  located  outside  of  O.  The  presenoe  of  the  two  spiky  estimates, 
however,  indicates  that  on  two  occasions, a zero  came  close  to  the  periphery 
of  O. 

One  of  the  major  drawbacks  of  this  technique  is  the  need  to  invert  a non- 
Toeplitz  matrix  (or  an  equivalent  operation)  in  order  to  evaluate  the  optimum 
filter  coefficients;  see  (140).  For  large  p,  this  is  a significant  numerical  prob- 
lem. We  therefore  attempted  to  convert  the  matrix  lSpm]  in  (140)  to  a Toeplitz 
matrix,  so  that  the  recursive  solution  in  appendix  A could  be  utilized.  We  first 
averaged  [Sf  m]  ? down  the  diagonals  and  left  the  right-hand  side  of  (140)  as  is; 
however,  we  lost  resolution  and  got  badly  biased  and  more  variable  spectral 
estimates.  Next  we  diagonally  averaged  [S/m]  P and  left  the  right-hand  side  of 
(140)  alone,  but  got  the  same  bad  effects.  Finally,  we  diagonally  averaged 
l Sfml  o ^ also  rePlaoed  the  terms  on  the  right-hand  side  of  (140)  by  the  ap- 
propriate averages,  but  again  to  no  avail.  Thus,  we  are  unable  to  significantly 
modify  (140)  without  dire  effects  on  the  spectral  estimate. 

Finally,  when  the  Burg  technique  is  considered  in  figure  10,  we  observe  the 
complete  absence  of  spiky  estimates;  this  is  due  mainly  to  the  guaranteed  loca- 
tion of  the  zeros  of  (A-10)  inside  O.  In  other  respects,  the  results  of  figures  9 
and  10  are  very  similar.  There  is  a small  bias  in  figure  10B,  with  the  peaks 
being  rounded  off  and  tho  valley  filled  in;  this  is  similar  to  figure  5 in  refer- 
ence 2. 

All  the  results  above  have  been  conducted  for  Gaussiai.  white  noise  jwjJ  in 
(198).  To  see  the  effect  of  the  statistics  of  |wj^[  upon  the  spectral  estimates, 
wo  changed  to  a uniform  distribution.  The  results  in  figure  11  are  virtually 
identical  to  those  in  figure  10.  Accordingly,  Gaussian  statistics  are  kept  for  the 
remainder  of  the  simulation. 
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6.2  BAD  DATA  POINTS 

By  virtue  of  the  results  of  the  preceding  subsection,  further  consideration 
is  limited  to  the  forward  and  backward  prediction  technique  and  the  Burg  tech- 
nique. The  first  example  we  consider  is  B ■ 4 bad  data  points  out  of  a total  of 
N = 40  data  points;  that  is,  in  each  of  the  100  realizations  of  40  data  pointB,  4 
points  (no  more,  no  less)  were  randomly  selected  as  being  bad,  and  the  corre- 
sponding values  of  Xfc  were  suppressed.  In  some  of  the  realizations,  the  four* 
data  points  may  have  been  oiose  together  (for  example,  10,.  12,  14,  ,16),  whereas 
in  other  realizations,  they  might  have  been  far  apart  (for  example,  1,  14,  27, 
40). 

The  resulting  speotral  estimates  are  given  in  figures  12  and  13.  . The  vari- 
ability in  the  skirts  is  less  for  th3  forward  and  backward  prediction  technique 
than  for  the  Burg  technique.  However,  the  spiky  nature  of  the  former  technique 
is  quite  evident  in  comparison  with  the  latter  technique.  Both  techniques  have 
suffered  a significant  loss  of  resolution  near  the  narrowband  peaks. 

The  reason  for  the  significant  degradation  in  performance  of  both  tech- 
niques is  that  although  only  B/N  = 4/40  (10%)  of  the  points  are  bad,  the  number 
of  valid  error  points,  N-p-Bp  in  (174)  and  (192),  can  decrease  significantly. 

For  example,  for  p = 4 and  spaced  bad  points  at  « 11,  M2  - 16,  M3  =»  21, 

M4  « 26  (see  (169) ),  we  have 

B = 20,  N ~p  - B - 16.  (200) 

p P 

On  the  other  hand,  for  contiguous  bad  points  at  =1,  M2  “ 2,  M3  ® 3,  M4  = 

4,  we  have 

B - 4 , N - p - B = 32.  (201) 

P P 

Thus,  anywhere  from  16  to  32  valid  error  points  can  be  achieved.  The  stability 
of  the  spectral  estimate  for  (200)  will  be  less  than  that  for  (201).  Generally, 
contiguous  bad  points  are  leBs  damaging  than  spaced  bad  points,  because  more 
valid  error  points  can  be  formed  when  the  bad  points  are  contiguous, 

One  of  the  points  of  the  above  example  is  that  4 bad  data  points  out  of  40  is 
rather  detrimental.  We  consider  now  N =»  100  data  points.  The  iirst  example 
of  interest  will  serve  as  a comparison  case  and  is  B » 0.  The  results  of  spec- 
tral estimation  for  the  forward  and  backward  prediction  technique  and  Burg 
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teohnique  are  given  in  figures  14  and  15.  respectively.  The  results  are  virtu- 
ally identical;  there  is  excellent  resolution  and  almost  no  bias  for  both  tech- 
niques . 

When  B is  increased  to  10,  the  results  in  figures  .16  and  17  are  obtained. 
Despite  10%  bad  points,  good  performance  in  terms  of  stability,  bias,  and  reso- 
lution is  attained.  The  number,  N-p-Bp,  of  valid  data  points  can  vary  from  46 
to  86;  however,  the  likelihood  of  realizing  as  few  aB  46  on  a random  basis  is 
very  remote.  The  Burg  technique  has  less-spiky  estimates  near  the  narrow- 
band peaks,  as  expected;  however,  it  is  more  variable  in  the  skirts  than  the 
forward  and  backware  prediction  technique. 


When  B is  increased  to  20,  the  results  in  figures  16  and  10  indicate  that  the 
Burg  technique  hue  more  variability,  but  is  less  spiky  and  has  better  resolution. 
The  same  conclusions  hold  true  for  B = 30  in  figures  20  and  21;  however, 
neither  teohnique  resolves  the  two  narrowband  peaks  for  thlB  many  bad  data 
points . 


s'  v.  [ ■ f i W'.i*  A ■ i,'  - 


Figure  4.  Yule-Walker; 


«Tii 


figure  5.  Yule-Walker,  Unbiased; 


Figure  6.  Least  Squares  of  Bos  and  J< 


um  Likelihood  of  Box  and  J< 
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Figure  11.  Burg,  Uniform  Noise;  N 
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Figure  18.  Forward  and  Backward  Prediction;  N = 100,  B = 20 
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7.  DISCUSSION  AND  CONCLUSIONS 

Several  methods  of  spectral  estimation  via  linear  predictive  techniques  have 
been  considered  for  a univariate  process,  both  with  and  without  bad  data  points; 
the  bad  points  can  be  regularly  spaoed,  randomly  spaced,  or  a combination. 

Two  particular  methods  have  been  found  to  have  better  performance  than  the  re- 
mainder, namely,  the  forward  and  backward  prediction  technique  and  the  Burg 
technique.  The  former  technique  tends  to  have  less  variability  on  the  skirtB, 
but  has  more  spiky  estimates  near  the  peaks  of  the  spectrum;  the  latter  tech- 
nique has  very  few  spiky  estimates.  Both  techniques  have  comparable  resolu- 
tion and  bias . 

Sinoe  the  best  ohoioe  of  filter  order,  p,  is  not  known  a priori,  it  1b  neoes- 
sary  in  practice  to  make  several  guesses  at  this  parameter  and  compute  some 
error  criterion  that  indicates  when  to  terminate  the  reoursion.  In  particular, 
Akaike's  Information  Criterion  (reference  22)  is  often  adopted  as  a termination 
procedure;  it  takes  the  form  (reference  1,  equations  (91)  and  (41)  or  reference 
22,  page  719) 

AIC  ■ in  (Relative  Error)  + |E-  (AIC  (p  = 0)  = 0)  , (202) 


where  Ne  is  the  "effective"  number  of  data  points,  and  is  taken  as  N-p  (or 
N-p-Bp  for  bad  points)  here,  at  the  p-th  stage.  The  value  of  p at  which  (202) 
is  a minimum  is  taken  as  the  best  estimate  of  this  parameter;  however,  criterion 
(202)  is  not  absolute,  and  the, user  can  actyuBt  It  to  fit  his  application  (reference  1, 
page  575).  A wide  range  of  values  of  p may  have  to  be  investigated  if  little  is 
known  about  the  true  spectrum  a priori;  an  upper  bound  on  p is  given  by  Akaike 
as  3N1/2  (Ibid). 

One  of  the  ramifioations  of  this  successive  guossing  procedure  is  that  for 
the  forward  and  backward  prediction  technique,  a different  pxp  matrix  [Snm]  P 
must  be  inverted  (or  an  equivalent  operation  oonduoted)  at  each  stage  (see  i 
(140)  and  (138) ) in  order  to  determine  the  filter  coefficients  and  minimum 
error,  (141).  Although  the  matrix  terms  can  be  updated  according  to 

X 4*  X 

S(P+D  _ N-p  (p)  p-H-m  p-KL-n  N-p+n  N-pi-m 
nm  ~ N-p-1  nm  2(N-p-l)  ’ d) 

in  addition  to  the  relations  in  (139),  the  size  of  the  matrix  tSnm]  ^ grows  with 
p,  and  the  solution  of  (140)  can  be  a time-consuming  procedure,  if  many  largo 
valuos  of  p mast  be  investigated.  This  fact,  coupled  with  the  fact  that  this 
spectral  estimation  technique  can  yield  spiky  estimates  and  on  unstable  recur- 
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sion  relation  (149),  leads  to  the  conclusion  that,  of  the  methods  considered,  the 
Burg  teohnique  is  the  recommended  prooedure  for  spectral  analysis  of  univari  - 
ate  processes.  A comparison  with  the  maximum  likelihood  teohnique  (reference 
23)  is  underway  and  will  be  documented  in  a future  report. 

The  solution  for  the  filter  coefficients  in  the  Burg  technique  is  accomplished 
recursively  as  shown  in  subsection  4.7  and  automatically  progresses  through 
successively  larger  values  of  p at  whioh  error  measures  (150t  and  (156)  are 
readily  calculated.  There  is,  of  oourse,  the  need  to  update  the  forward  and 
backward  residuals  via  (153),  and  the  calculation  of  oross-galn  gp  in  (155),  both 
of  which  t.-ke  time  to  effeot.  But  the  effort  required  actually  decreases  as  p 
increases,  since  fewer  terms  are  involved  in  (153)  and  (155);  in  exchange,  tho 
stability  of  the  estimates  also  decreases. 

FORTRAN  programs  for  the  Burg  teohnique,  both  with  and  without  bad  data 
points,  are  given  in  appendix  J.  Some  representative  execution  timcB  on  the 
Univao  1108  for  the  computation  of  the  filter  coefficients  (SUBROUTINE  BURG) 
are  given  in  table  3,  where  N is  the  number  of  data  points  and  PMAX  is  the 
maximum  order  of  filter  considered.  The  times  are  approximately  linearly 
proportional  to  N and  PMAX.  The  execution  time  for  the  evaluation  of  the 
power  density  estimate  itself  is  governed  by  the  FFT  teohnique  employed  to 
evaluate  (167)  (SUBROUTINE  POWERS). 


Table  3.  Execution  Times;  No  Bad  Data  Points 


N 

PMAX 

Time  (sec) 

100 

10 

0.038 

100 

20 

0.073 

1000 

10 

0.33 

1000 

50 

1.78 

10000 

GO 

17.9 

10000 

150 

48.4 
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The  presence  of  bad  data  points  Is  easily  aooommodated  In  the  Burg  tech- 
nique, as  shown  In  subseotion  4.7.  If  the  bad  data  points  are  contiguous,  the 
loss  in  stability  of  the  estimates  Is  not  as  great  as  when  the  bad  data  points  are 
spaoed.  The  worst  possible  locations  of  bad  data  points  oocur  when  the  olosest 
spaoing  is  >p  + 1,  sinoe  each  bad  data  point  causes  the  loss  of  p + 1 valid  error 
points.  Interpolation  of  spaoed  bad  data  points  has  proven  poorer  than  the  tech- 
nique utilized  here  (of  ignoring  bad  points)  when  the  speotral  oontent  of  the  input 
process  extends  fairly  dose  to  the  Nyquist  frequency  (2  A)-^.  Sinoe  the  exact 
extent  of  the  Input  speotrum  is  unknown  a priori,  interpolation  can  be  a damag- 
ing procedure  in  some  cases. 

The  spectral  estimation  technique  investigated  here  Is  particularly  advan- 
tageous for  short  data  segments,  where  other  methods  are  inapplicable.  For 
example,  if  a piece  of  equipment  fails  frequently,  short  disjointed  pieces  of  data 
may  be  all  that  are  available.  Or  if  a prooess  is  nonstationary,  it  may  be  nec- 
essary to  out  the  total  data  record  into  small  segments  in  each  of  which  it  is 
believed  that  conditions  are  substantially  stationary.  For  longer  data  records, 
where  standard  FFT  techniques  can  be  applied,  it  has  been  recommended  that 
both  speotral  estimation  procedures  be  applied  and  the  results  plotted  together 
to  glean  maximum  information  about  the  true  speotrum  (Bee  referenoe  12) . 

This  seems  particularly  useful  when  some  pure  tones  are  present  in  the  Input 
data;  the  standard  FFT  technique  is  ideally  suited  for  the  analysis  of  pure  tones 
or  very  narrowband  components. 
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Appendix  A 

RECURSIVE  SOLUTION 


If  we  employ  (52)  in  (46),  there  results 


Now  define 


P 

E 

k=  0 


R 


jf-k 


4,  , 0<  H £ p. 

lo 


pr  59  ■ . 0 ^ k < p , 


(A-l) 


(A-2) 


where  the  dependence  of  the  ooeffioients  on  the  order  p in  (31)  is  indicated  ex- 
plicitly. Then  (A-l)  becomes 


- m 

* 

»•  - 

R R , ...  R 

1 

l/oW 

o -1  -p 

oo 

R,  R 

-a?> 

0 

1 o 

1 

• 

a 

. 

• , 

R R 

• 

-a<» 

0 

P o 

4 *• 

L p J 

- - 

(A-3) 


where  the  matrix  R is  Hermitian  and  where  we  have  also  indicated  that  the  real 
quantity  c00  is  dependent  on  p;  see  (47)  and  (51).  Equation  (A-3)  constitutes 
p+1  linear  equations  in  the  p+1  unknowns  a , ....  ajlp) , l/o^P)  . 


The  solution  to  (A-3)  oan  be  obtained  recursively  as  follows  (see , for  ex- 
ample, reference  11  or  reference  24,  appendix  B): 


= R /R 
1 o 


R 

o 


- R a(1) 
-1  al 


(A-4) 
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for  p > 2 : 
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c(P"1) , (A- 5) 
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, k - - 1 , 


p-k 


(A- 6) 
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o *?  k o p P k k „(P-1)V  P ) 
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k»  i 


k"3 1 


oo 


(A-7) 


The  last  step  in  (A-7)  is  obtuined  by  substituting  (A-6)  and  employing  (A-5).  It 
is  very  important  to  notice  from  (A-tt)  that  onue  9^  is  specified,  all  the  p-th 
order  filter  coefficients  can  be  calculated  from  (p-l)th  order  coefficients . The 
same  is  true  of  (A-7) . 


If  we  uBe  (A-2)  and  (53),  the  maximum  entropy  spootrum  in  (55)  oan  be  ex- 
pressed as 


G (f) 
o 


4/0®> 

OO 


V -to) 


V ,£,<2T‘ 


a£'  exp(-i2»rfkA) 

k=0 


(A- 8) 


The  similarity  in  form  to  (14)  will  be  complete  wher.  it  is  shown  (in  (67) ) that 
1/oto)  is  the  minimum  value  of  the  average  magnitude- squared  error  for  a p-th 
order  predictive  filter;  therefore  must  be  positive  for  all  p„  for  non-negativa 
definite  R.  Equation  (A-7)  offers  a recursive  calculation  of  the  average  error;  It 

oan  be  started  with  R0  • (In  fact,  (A-5)  through  (A-7)  can  be  used  for  p > 1 

when  that  starting  value  is  used. ) 


A-2 


1 
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Since  must  be  positive  for  all  p,  (A-7)  indicates  that 

< 1 for  k ■»  1,  2,  p.  (A-9) 

This  is  equivalent  to  having  all  the  zeros  of 

V a?>  a-k  , <A-10> 

k“o  k 

(where  the  remaining  coefficients  are  determined  via  (A-8) ) inside  the  unit  cir- 
cle, G,  in  the  complex  z -plane;  see  reference  1,  page  667,  Therefore 
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Appendix  B 

EVALUATION  OF  MAXIMUM  ENTROPY 


The  optimum  spectrum  is  given  by  (33)  and  (37).  The  maximum  entropy 
then  follows  from  (30;  as 


/ 

^l/A 


l 


Ent  ■?.  4|  df in  G^(f)  = -A  I d£tjCn*y(t‘)  +j2aV  (01  2 ^ + £2  . 


Consider 


^1  •"  - A f dfj2n  | ^ exp (i2  rr  f k A)  j . 

A/a  (k=»  o ) 


(B-l) 


(B-2) 


Lotting  /,  « uxp(i2  ir f A)  and  using  (38),  (B-2)  becomes 

= - st/t  ^ • <B-3) 

where  y denotes  counter olookwise  Integration  around  me  unit  oircle,  O,  in  the 
oomplex  z -plane.  Now 


B(z)  » 
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k=  0 
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where,  from  (A-14),  zero  locations  |o^|  satisfy 

|ok|  > 1,  all  k; 


(B-4) 


(13-5) 


that  is,  all  the  zeros  of  B(z)  lie  outside  O.  (There  can  be  multlplo-oroer 
zeros  in  (B-4). ) Also  assume  p for  now.  Then  (B-3)  can  be  expressed  as 


B-l 
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But 


" ' ST  /-T  p"  “p  + v] 

L k«l  J 

= - Fine  t (»-  oj  • 

L k=l  J 


(B-6) 


in  (z  - o^)  **  in  (-  0^  + in  (1  - ~-) 
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(~\  - ...  for 

z 

\\) 

\ 

< 1; 


(3-7) 


that  Is,  expansion  (B-7)  converges  for  Izl  < |ok|, . But  since  lo^l  >1,  the  region 
of  integration  in  (B-6)  remains  in  the  convergence  region  of  (B-7).  Therefore, 
the  integral  in  (B-6)  is 


Then  from  (B-6)  and  (B-4) 


= - pn  «p  + £ (*  °k)  " ~ pp  TI  ("  °k)j 


* - in  B(0)  = - in  aQ  . 


(B-9) 


And  from  (B-l)  and  (B-4) 
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Now 
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in  ('  0 ' ^ ' fe)  ' 


...  for 


°kz 
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(B-ll) 


that  la,  expansion  (B-ll)  converges  for  Is)  >^f.  But  since  lokl  >1,  the 
region  of  Integration  In  (B-10)  remains  in  the  convergence  region  of  (B-ll). 
Therefore,  the  integral  in  (B-10)  is 

±#**e  ■ <)■  ifih  (••»-#•(#)'■  -I  ■ " ('5!: 


(B-12) 


Then  from  (B-10)  and  (B-4) 


- [*"  “p  * JB "5?]  - - ii^)] 


- in  B* (0)  = ••  in  a< 


(B-13) 
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Combining  (B-9)  and  (B-13)  in  (B-l),  there  follows  for  the  maximum  entropy 

Ent  = - jfi  n f °<0 1 2 " in  > (B-14) 

where  we  have  also  employed  (52).  (For  p = 0,  a separate  derivation  yields 
(B-14)  also.)  Recall  from  (51)  that  o00  is  the  upper-left  corner  element  of  R*1, 
where  R is  defined  by  (47) . 
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Appendix  C 

IMPLICATIONS  OF  ASSUMPTION  OF  WHITE  SPECTRUM 
FOR  MINIMUM  ERROR;  KNOWN  CORRELATION 


We  define  the  crossoorrelation  function  between  minimum  error  7 and  Input 
x in  figure  1 as 


(C-l) 


Substituting  (59)  and  utilizing  (1) , this  becomes 


C,  - £ Z R,  , alU 
1 n f-n 

n»  0 


(C-2) 


Now  from  (64)  and  (66),  we  can  express 


R ■ * - — — t . 
0 

00 


(C-3) 


Thus,  (C-2)  immediately  yields 


Cf  " 


-l/o  , j2  a 0 
oo 


0,  1 < JL  < p 


(C-4) 


that  is , minimum  error  value  7^  is  unoorrelated  with  the  past  p inputs  x^-i, 

• • • i ^k-p ■ 

Now  using  (59)  and  (C-l),  the 'autocorrelation  funotion  of  the  minimum 
error  is 

P ____  P 

E.  s 7.  7*  = T.  H*  r.  x*  t - V 1*  C.  . all  J2  . (C-6) 

t k k-i  *-*  n k k-f-n  n f+ n 

n=*  0 n*0 


r iV, , ..1  I .'Vtr4  '*'*t  ' 
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In  particular,  using  (C-4), 


But  from  (C-2)  and  (66), 


E,  - **•  C . 

1 p p+1 


C . - T.  It  R . , ■ - R 

p+1  ^ n p+l-n  ] 

n-0 


P o 

. . , r — b . , . 

P+1  „.l  °oo  ptl-" 


Therefore,  assuming  E^  ■ 0 is  equivalent  to  assuming  Cp+i  « 0,  (that  is,  mini- 
mum error  7^  unoorrelated  with  input  x^-p-i),  which  in  turn  is  equivalent  to 
requiring 


r » ~ y 2 r ■ r t h . <c-8) 

p+1  ^ 0 p+l-n  **  n p+l-n  ' ' 

n=l  oo  * n-1  r 

This  relation,  whioh  may  not  be  true  for  the  actual  process  |xjc[ , is  a dlreot 
result  of  assumption  (70);  the  quantity  Rp+i  in  (C-8)  is  really  an  approximation 
to  the  true  (unknown)  correlation  value. 


Next  from  (C-S), 


E9  ■ X 1 + X 4.9  • 

i P“1  p+1  P P+2 


Assuming  Eg  “ 0 (in  addition  to  Ej  ■ 0)  is  equivalent  to  also  assuming  Cp+2  « 0, 
whioh  in  turn  from  (C-2)  and  (66)  requires  that  we  approximate  aooording  to 

V'^r  V-n  ■ k *»  V-n  • <C-10) 

n«l  oo  n-1 

Continuing  in  this  way,  it  follows  that  assuming  white  noise  for  {7^} , that 
is,  assuming 

E,  - 0 for  j >1,  (C-ll) 

is  equivalent  to  assuming  that  Ct  - 0 for  SL  > p + 1;  that  is,  the  minimum  error 
is  unoorrelated  with  all  past  inputs.  There  follows  the  approximations 


*v  : • ‘ , \ . -p,  V*..  UW.  .V!  L-~  * : 
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«i**  E E *nR,-»IOr 

n*l  00  n«l 


(C-12) 


This  recursion  relation  (starting  with  known  values  R^,  R2.  ....  Rp)  oan  be 
considered  to  be  an  extrapolation  of  the  known  correlation,  values  into  regions 
where  they  are  unknown. 


If  we  augment  (C-12)  according  to 


R ■ R*  forj2>  p+1  , 
then  it  can  be  shown  that  the  spectrum  defined  by 


(C— 13) 


A £ R,  exp(-12irftfA)  = 
|o-« 


A/o 


oo 


1 - ]j£  a exp(-i2irfnA) 
n«l 


(C-14) 


whioh  is  identical  to  (71).  The  transform  in  (C-14)  converges  if  | Rp  | deoays 
with  increasing  (i  \ , that  is,  if  B(z)  of  (56)  hac  no  zeron  inside  0 . 


C-3/C-4 
REVERSE  BLANK 
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Appendix  D 


STABILITY  OF  RECURSION  RELATION 


The  reoursion  relation  for  approximated  correlation  valuer  ■ f is  given  in 
(C -12)  and  (75)  as 


% - E \ R^k  for  i > P + 1 
k«l 


Therefore, 


CD  | * 

u««  E R,  z"  - E V"*  2 Vk 

f«p  + l k-1  /“P  + l 


z'i+k.  (D-2) 


<»  Ij-t,  u>  , H 

E R,v  - E *,»  - E »j 

i»p  + l j»p+l-k  j«p  + l-k 


«>  , 

+ £ J * Vk(a)  + U(a> , (D-3) 

)«p  + l 

where 

Vk(z)  * Rp+l-k  z_(P+1’k)  + Vk-1(z) , k £2 ; ^(a)  - Rp  . (D-4) 

Vk(a)  involves  the  starting  values  Rp+i_k>  ....  Rp  for  1<  k<  p,  Employment 
of  (D-3)  in  (D-2)  yields 

P . P . 

U(a)  “ E \ z Vfc(a)  + u(z>  E \ . (D-5) 

k®  1 k«l 


I,,,,  , 
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or 


U(z) 


h \ *'k  vk<z)  £ cko vk« 

k»l  _ k-1 


(D-8) 


E V 

k=»  1 


-k 


E o 

k=0 


ko 


I 


where  we  havo  utilized  (66).  In  order  that  recursion  p-1)  be  stable,  the  de- 
nominator of  (D-6)  must  possess  all  its  zeros  within  the  unit  circle  O in  the 
complex  z-plano.  Therefore,  B(z)  of  (56)  must  possess  all  its  zeros  outside 
O if  recursion  p-1)  is  to  be  stable.  This  is  guaranteed  by  the  results  in 
(A -9)  et  seq. 


D-2 


‘ , l 


A'  ***  - i.1 

■ k-v  rtk.  i* 
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Appendix  E 


IMPLICATIONS  OF  ASSUMPTION  OF  WHITE  SPECTRUM: 
UNKNOWN  CORRELATION 


The  minimum  error  sequence  is  given  by  (96)  and  (101)  aB 


T.  a H * x.  , all  k, 
k ^ n k-n 
n«0 


(E-l) 


The  sample  autocorrelation  of  {TjJ  is  defined  here  aa 


f.  r.  « y\  k s, . (e-2) 

* N k k-P  **  m n /+n-m  ' 

k m,n=>  0 

using  (E-l)  and  (98).  The  sample  Bpeotrum  of  { is  defined  here  as 

Hr(f)aA  f;  F#  exp(-i2irfiA)  - Hx(f)|A(f)|2,  Ifl  (E-3) 


jfa  -uu 


where  we  have  employed  (E-2)  and  (107)  and  defined  the  sample  spectrum  of 

ixkl  aB 


Hx(f)2A  E S|  exp(-i2nfiA),  If l<~- 

/ ■>  -to 


(E-4) 


Therefore,  (E-3)  yields 


H (f) 
x'  ' 


Hf(f> 


A(f)| 


2-  1,1  <U' 


(E-5) 


Now  we  will  assume  that  the  sample  speotrum  of  {T^}  is  white;  that  is,  we 
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l 4W* 


HT(f)  s Kd , Ifl  < ~, 


(E-6) 


where  K is  a constant.  We  then  adopt  an  estimate  of  the  sample  spectrum  of 
sequence  |xn| ^ according  to 


H (f)  3 
x' ' 


H?(f) 


Kd 


A(f)  • 


f exp(-12?rfnd) 

n™  0 


•,  Ifl  < 


2d  ’ 


(E-7) 


and  adopt  a scaled  version  of  this  quantity  as  a spectral  estimate  of  process 

lxni  : 


G (f)  9 
x' ' 


^ ft  exp(-12rrfnd) 
n»0 


(E-8) 


The  white  assumption  in  (E-6)  forces  ub  to  assume  that 

F(  ■ 0 for  £ 0 , 


(E-9) 


as  (E-3)  shows.  In  order  to  see  what  this  implies,  we  utilize  the  definition  of 
the  sample  crossoorrelation  in  (109),  along  with  (96)  and  (98),  to  obtain 


v,  - 1 

k n=  0 


(E-10) 


Use  of  (101)  then  shows  that 


Dj  “ 0 for  1 < < p 


(E— 11) 


Meanwhile,  the  sample  autocorrelation  in  (E-2)  can  be  written  in  the  form 


pf  " £ D/+n’  alli' 

n-  0 


(E-12) 


E-2 


»lf4 

»lil  - # - ii -l.'1- 
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t 


i 


upon  employment  of  (E-10).  There  Immediately  follows  from  (E-9),  (E-ll), 
and  (E-12) 


F = V D . , » 0 . 
1 p p+l 

(E-13) 

But  then  (E-9)  and  (E-10)  indicate  that 

P 

s = V S s ^ , 

p+l  *-»  n p+l-n 

n=  1 

(E-14) 

•> 

where  |a^|  £ are  the  solutions  of  (102).  But  relation  (E-14)  may  not  be  true 
for  the  quantity  Sp+i  actually  obtained  from  data  jxn|  ^ via  (98).  Thus,  as- 
sumption F'i  » 0 is  forcing  us  to  assume  that  Sp+^  can  be  obtained  via  (E-14) 
and  (102),  when  {s^  }P^  are  obtained  from  (98). 

Next  from  (E-12)  and  (E-ll), 

F„  » T . D ■+  Tt*  D 0 . 

2 p-1  p+l  p p+2 

(E-15) 

Assuming  F2  a 0 (in  addition  to  Fj  ® 0)  is  equivalent  to  also  assuming  Dp+2  =>  0, 
which  in  turn  from  (E-10)  requires  that 

P 

S O v ft  s 
p+2  n p+2-n 

Tl~  1 

(E-16) 

Continuing  in  this  way,  it  follows  thut  assuming 

F(  = 0 for  JL  >1  (E-17) 

is  equivalent  to  assuming  Df  = 0 for  i £p  + 1;  that  is,  the  minimum-error  se- 
quence is  uncorrelated  (on  a single  member  function  basis)  with  all  past  inputs. 
There  follows  the  estimates 


P 

S(  = L S|_n'  i>P  + i.  (E-18) 

n=>  1 

Stability  is  disouBsad  in  (111)  et  seq. 


E-3/E-4 
REVERSE  BLANK 
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Appendix  F 

BOUND  ON  CROSS-GAIN 


The  value  of  the  orosB-galn  gp  in  (155)  can  be  written  as 


f(p~l)  b(p-l)* 


The  first  factor  in  (F-l)  is  of  the  form  of  a correlation  coefficient  of  the  (p-1)- 
th  order  forward  and  backward  sequences  and  can  never  exoeed  unity  in  magni- 
tude (by  Schwarz's  inequality).  The  seoond  factor  in  (F-l)  is  almost  always 
very  close  to  1:  let  the  pair  of  sums 


{ ^ 

Vn«p+1 


P(P-D| 


and 


N 

£ 

=*p+i 


Vi 


' | e |a  and  A(l+r)|  , 


(F-2) 


where  £ >0  without  loss  of  generality.  The  seoond  faotor  in  (F-l)  then  equals 
(l+r)' 


l/2 


T+^r 7T~ * whio^  ls  never  than  1 and  ts  tabulated  below.  Thus,  gp  in 

(F-l)  is  virtually  Identical  to  the  correlation  coefficient  of  the  forward  and  back- 
ward sequences,  sinoe  r is  near  zero  with  high  probability. 


Table  F.l  Second  Factor  in  (F-l). 


r 


0 


.1  .2  .3 


.4  .5 
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i : i: 


Appendix  G > 

\ \ ; , . ' 1 

CLOSENESS  OF  ERROR  MEASURES 


Two  possible  error  measures  for,  the  Burg  technique  were  presented  in 
(150)  and  (156).  For  p ■ 0,  employing  (154)  and  (152), 


r<0)  ■ i £ k* 

n=  1 


I •,  : 


Comparing  this  result  with  (151),  we  find 


F<°>  . p<°> 


Thus,  the  two  error  measures  are  identical  for  p = 0. 
Next  from  (160)  and  (161) 


p«  = P<°> 


(i-Kl,f)-(i-K1)r)s  2|>j2 

' ' ' n=l 


whereaB  from  (156),  (160),  (155),  and  (152), 


Den(l) 

2(N-1) 


- k - k«r) 

- (- -KT)^i>  s (l^l2  -tef) 


(G-l) 


(G^2) 


(G-3) 


kw0 

l 

± 

X1 

2 

+ X2 

2 + 

|X3 

2 + ...  *| 

i 

XN-2 

2 

+ 

Vl 

1 4 

2 

XN  ‘ 

^ N- 1 

G*1 
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But  now  reference  to  <151)  and  (G-3)  reveals  that,  for  N-l  large. 


(G-b 


Continuing  with  (150), 


(G-e 


And  (156),  (160),  and  (156)  combine  to  yield 


'?  - (‘  - I4T)  SSSH  - Is*  f ) ohr  £ (Iff  * ICO 

n“3  (G- 


But  from  (154), 


-r-sferW-ICl') 

n-2 


(G-8 


Comparing  (G-7)  and  (G-8),  we  see  that,  for  N-2  large, 


F«» 


(G-9 


(G-lo 


(-1st) 

Then  employing  (G-5)  and  (G-6),  we  have 

P«2,»(l-|aff)p<1)-P(2>, 

which  is  the  desired  relationship.  In  general,  for  no  bad  data  points,  we  have 

F<P)  * pk)  for  N-p  large  . (G-i: 

Numerical  computations  have  borne  this  result  out,  with  the  two  quantities  not 
having  any  ordered  relationship;  that  is,  either  quantity  oan  be  larger  (or 
smaller)  at  different  stages,  p.  (G-9)  generalizes  to 

* ^1  - |a^|2^  F^”l)  for  N-p  large.  (G-1S 

G-2 


C - *K.  ■ -.4 ' '-pH  m- 


Tft  6303 


Appendix  H 

SCALE  FACTORS  IN  SPECTRAL  ESTIMATES 


Instead  of  using  a unity  value  for  the  average  minimum  error  or  residual 
power  In  the  numberator  of  (167),  we  oould  use  the  value  given  by  (156),  Then 
our  speotral  estimate  would  be 


8x® 


A F 

- £ 4P) 

k»l 


(P> 


exp(-i2rr  fkA) 


Ifl  < 


1 

2A 


(H-l) 


An  alternative  approach  is  to  use  an  arbitrary  Boale  factor  K and  ohoose  it 
so  that  the  area  under  the  speotral  estimate  is  equal  to  the  sample  power  (l&l), 
as  suggested  under  (108);  that  is,  set 


Sx® 


AK 


1 - £ a^exp(-i2fffkA) 
k°*  1 


- If  I < -a- 
2*  111  2A* 


and  foroe 


i 

2A 


“ N 

f df6,«-p«-A  £ jx„|2. 

J n»l 

~2A 


Subotituting  (H-2)  in  (H-3),  and  using  (159),  we  havo 


,(0) 


1 

2A 

■/. 


df 


AK 


K 


2A 


(H-2) 


(H-3) 


p 

2 P , 

, , 2 

1 - 2 \ exp(-i2fffkA) 

n 1 

1 - a(m) 
I m 

k-1 

,m»l  1 

(H— 4) 


\ 

'J 

1 

i 

* 


H-l 


I 
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The  last  step  in  (H-4)  is  proven  as  follows:  from  (A-8)  and  (29),  we  know  that 


1 

2A 


J 


1 • df  A'.TVf  -ta  Jr;  , 


P 


■ I . 1 1 1 M - , - i ■ I I:  :■> 


-JZ  1 - £ exp(-i2fffkA) 

k«l  I-  ' 


i . i . i ; 


0 00 


(H-5) 


■I.  ■,!  I i-  l I 

But  from  (A-7) , 


. • • ! i 


ill'll  ■■  ;t|  ■ ' i : i>  ! li.'  .)  if 

I 1 1 ■ 1 , ;l  : , I ■ , . I 

■ . i ; 1 1 i i * 1 ■ i i i : i 


R O4’-1*  "!• 


R 0«  . ° 


O 00 


1 - a 


(P)|2 


m-1 


tj  H - *' 


(H-6) 


* I 


where,  we  have  employed  R0  ofa  ■ 1, , The  relationship  in.  (H-4)  holds  when  the 
filter  .coefficients  are  determined  via  (148).  ,, , 


: l • 


Therefore,  (H-4)  yields,  with  the  aid  of  (150). 


.«* 


K » Pi.  n 1 - K 


| m-1 

and  the  estimate  (H-2)  becomes 


m 


M 


(H-7) 


! I 


6x(f) 


AP 


<P> 


i.  - Y)  ^ exp(i-i2ff  fkA) 
k»l* 


(H-8) 


The  very  close  similarity  of  values  between  the  alternatives  (H-l)  and  (H-8)  is 
made  evident  by  the  results  Of  appendix  ti,  lh  particular  (fr-11),  Thus,  there  is 
virtually  no  difference  between  estimates  (H-l)  and  (H-8),  for  no  bad  data 
points . 


i 


H-2 


"i>  ! I1  . 'I 

Appendix  I 


, -:l.  1 ' 


I I 


. > <; r 


, u.  .BIASEpNTESS  pF^RG'S, CORRELATION  ESTIMATE 

1 1 i ■■'! ! I '!■:  1»  <i  i i 1 1. 


For  the  Burg  technique  with  p = 1,  N ■ 3,  we  find  from  (162)  and  (144)  that  i 
(for  real  data)  . -i  . ' f . ' u 


~ 2 xa  (xi + xa)(x?  + x2  + xs) 

R o — * ■ — !■■■'■ • t I 

n n a n 1 • 


1 3 


2 2 2 

V2V*3 


(M) 


\ n 


The  mean  of  this  random  variable  depends  on  more  than  just  X2X-L  (*  x3x2 ) ; in 
fact,  it  depends  on  the  third-order  joint  density  of  (Xp  x2,  x3).  As  an  exam- 
ple, let 


1 " X2  a “ v?  (u  + v)»  xg  a v» 


(1-2) 


where  u and  v are  independent,  zero-mean,  unit- variance,  Gaussian  random 
variables.  Then  Xp  x2,  x3  are  zero-mean,  unit-variance,  Gaussian  random 
variables  with 


X2  X1  51  X3  X2  * — ; *3  *1  B ° ’ 


(1-3) 


Employing  (1-2)  in  (1-1) , we  obtain 


•5  1 1 (u  + v)2(3u2+  2uv  + 3v2) 

Ri  " ij  ? - * ’ 


2x  . "2" 
u + uv  + v 


(1-4) 


Therefore, 


£ 11  1 f f j , ( u2  + v2  \ (u  + v)2  (3u2  + 2uv  + 

H1  ■ is  6 27  J!  du  dv  ex»\-— ~ J t—  ■ 2 

Jym  ' u + uv  + v 


.2  2 


3v2) 


TCO 


H 


V , 


MtOwrlrtWMtl-AH, 


1 
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+ii  i r 
■V2  8 2ff  J 


I dr  r3  exp(-r2/2)  f d,  M 

l J cz  + cs  + sz 

u -rr 


where  we  have  changed  to  polar  coordinates  and  let  C b cob  6,  S b sin  6.  The 
Integral  on  r In  (1-6)  Is  2,  and  the  Integral  on  6 Is  4irf2  - IA  . 


Therefore, 


<-94e4)- 


whloh  Is  not  equal  to 


x2  X1  " *3  X2  “ * jz  ’ 


''wfi 


r 
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Appendix  J 

FORTRAN  PROGRAMS 


The  programs  in  this  appendix  are  written  for  real  data,  but  may  be  readily 
generalized  to  complex  data  by  means  of  the  general  equations  in  the  main  text. 
From  (H-8)  and  (H-7) , for  real  data,  the  speotral  estimate  ia  given  by 


8 (f)  - 
x' ' 


* 77  1 " \ 
k-1  ' K 


(k)' 


,(0) 


1 - a^  exp(-i2  rrfkA) 
k“l 


2*  ,fl<*’ 


(J-l) 


Let  frequency  increment 


JL  1 JL  fN 

Af  S JA  " J72  2A  “ 572  ’ 


(J-2) 


where  fN  is  the  Nyquist  frequency,  and  J is  an  integer.  Then,  using  (H-3)  and 
the  real  behavior  of  the  data. 


_1_ 

2b 


,(0) 


. 


df  G (f)  a 2np<0)  [J 
* k-1 


J/2 

Z)  • 

ma  0 


■m 


2p(0) 

J 


1 - 

2 •£' 

1 k®  1 

p 

r ...2 

. J/2 

77 

1 - 

v ! 

E T 

k»l 

: ] 

' m = 0 

‘m 


1 - £2  exp{-i2  ffmk/J) 
k“  1 


J/2 

= T2  « p . 

mm 


2 ' *-*  mm 
m®  0 


(J— 3) 
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where  |«m|  is  a set  of  Integration  weights  (for  example,  trapezoidal).  So  we 
oan  compute  (independent  of  time  inorement  A)  the  quantity 


2 ’ 
j n 

k-l 


- £ •?  ’ 

k-l 


exp(-i2jrmk/J) 


for  0 < m < 


J 

2 * 
(J-4) 


which  represents  the  fraotional  power  In  the  frequency  band 


that  is, 


/m  - a 

m+2 

\ JA  ’ 

JA 

J/2 

P 

T'  i 

m 

m-0 

-JEi 

p(0) 

(J-5) 


(J-6) 


if  estimate  Gx(f)  in  (J-l)  has  been  sampled  finely  enough  (that  Is,  large  J).  The 
denominator  of  (J-4)  is  reoognized  as  a J-point  FFT  of  p+1  nonzero  numbers; 
henoe,  J should  be  chosen  as  a power  of  2 for  speed  purposes.  The  programs 
below  yield  the  fraction  of  power  in  frequency  bands  of  width  (JA)-1,  if  J is  an 
Integer  large  enough  that  the  spectral  estimate  (167)  or  (H-8)  is  adequately  sam- 
pled to  keep  track  of  itB  peaks . 


NO  BAD  DATA  POINTS  (SUBSECTION  4.  7) 

The  data  generation  is  accomplished  via  funotlon  IRAND,  which  generates 
integers  uniformly  distributed  over  (0,  233-  i)>  by  RAND,  which  generates  num- 
bers uniformly  distributed  over  (0,1);  and  by  TINORM,  which  generates  zero- 
mean  unit-varianoe  GausBian  variables.  The  FFT  used  below  is  that  presented 
in  reference  25. 
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t SPtCTHAL  ESTIMATION  USER!  CHANBE  LJNt  13  ANO  REPLACE  LlNtS  17-31 
C N a NUMBER  OF  DATA  POINTS 
C XU)»t«"X(N>  B INPUT  OATA 
C PM AX  B MAXIMUM  ORDER  OF  FILTER 
C PBtiT  s BEST  ORDER  OF  FILTtR 

C AU)»...»A(PBEST)  B PREDICTIVE  FILTER  COEFFICIENTS 
C PROD  ■ PRODUCT <1-A(P>**2)  FOR  PH  TO  PBCST 
C HHu(1)m..»RH0(PmAX)  B NORMALIZED  CORRELATION  COEFFICIENTS 
C J * SIZE  OF  FFT  (MUST  BE  A POWER  OF  «) 

C XXU>»,,..XX(J/a+l>  b FRACTIONAL  POMEHSi  FROM  DC  TO  NVfiUIST  FREQUENCY 
C C0(1 ) i • • • rC0(U/4+l)  a QUARTER  COSINE  TABLE 
C Y AND  YY  ARE  REQUIRED  AUXILIARY  ARRAYS 

PARAMETER  Nb  100*  PMAXaiO»  ja204B»  J4&BU/4+1 

integer  pbe&t 

DIMENSION  X(N) » T (N) » A (PMAX J |RHO (PMAX) iXX(O) r YY ( J) «C0(U41 ) 

L INPUT  DATA  IN  X(l> # . ..#X(N) 

DkF INE  IRANOB 148**18+ ( ( i-SION (1 * I*»***B) ) /«) *34359730367 
Dth INE  R aNObFlO AT ( 1 ) t 34389738367 . 

I.', 6281 

NSTAMTsN+400  « A’XLL  DISCARD  INITIAL  400  POINTS 
XX(1)B0, 

XXU)BO, 

XX13)B0, 

XX(4)B0, 

DU  11  LaBfNsTART 
lalRANO 

XX (L  > b2, 7607«XX (L-l ) -3 . »106*XX ( L-2) +2.6b3b»XX (L-3 1 - 
*0 , 923B*XX (L-4 ) *T  XNORM (R*ND . $11 ) 
n continue 

UU  12  UliN 

IS  XUJBXAII+NSTART-N) 

PRINT  li 

l format (/*  input  uatao 

PRINT  4,  (X(I)'ial'N) 

c evaluate  PhEuxcTive  falter  coefficients 
Call  BUR4(N»PMAX»X*YiPl9eST,A»PR0D,RH0i 
Print  v.  x(n ) 

9 FORMAT (/ f MEAN  s'fLl4,ti) 

PRINT  10»  T < N ) 

10  FOHMAT ( * STANDARD  DEVIATION  a'*EI3.B) 
print  2.  PmeSt 

d FOKMAT(/J  Petsr  at*l3) 

print  a, 

o formats  predictive  filter  coefficients!') 

PRINT  4,  (A(l)iiai.pBEST) 

4 FoKMAT(3E20,B) 

PRINT  b*  PROD 

b FuHMAT (/ » PRODUCT (1-A(P)«*2)  :'«E13.6) 

PRINT  6. 

t FORMAT*/'  NORMALIZED  CORRELATION  COEFFICIENTS! • ) 

print  4,  irho<i>*i«i*pmax> 

CmuL  QTNCOS(COiJ) 

c evaluate  fractional  powers 

Call  P0W£R5(PBEST»AiPRQUfU»XX#YY»C0) 

PRINT  7, 

7 FORMAT (/*  FRACTIONAL  POWERS!') 

L 30/241 

PRINT  b.  (XX(X) »i*l*L) 
o Format (2Xr loEl3«6) 

END 


> *Zih  Vi*4'i4*V.  di* 


orr  r.  -j  ro 


subroutine  burb(n#pmax#x*y#pbest#a#prod#rho>  e t feb  int 
c this  subroutine  computes  the  predictive  filter  coefficients 

C N a NUMBER  OF  DATA  POINTS  I INTCSER  INPUT 
C PMAX  ■ MAXIMUM  ORDER  OF  FILTER  I XNTESCR  INPUT 

c xu?»x(«)*,,.»x<n>  a data  aRRAv  on  input f altered  on  output 

C ON  OUTPUT#  XU)’X(8)#>»>#X(PMAX)  a AUlFMAX)#A(8)PNAX)»...,A|PMAX)PNAX) 
C Yll) # Y (8) # «« • # Y(N)  > AUXILIARY  ARRAY | SCRATCH  INPUT 

C ON  OUTPUT#  y (I) • V (2) > • • t • Y(PMAX)  a #A(8)8> A (PNAXIPMAX) 

C ON  OUTPUT#  X(N>  a MEAN#  ANO  YIN)  a STANDARD  DEVIATION  OF  INPUT  DATA 
C PBEST  b BEST  OROER  OF  FILTER)  INTEBER  OUTPUT 
C A< 1 ) » A(8) » , , . # A (PBEST)  a PREDICTIVE  FILTER  COEFFICIENT  ARRAY  a 
C All) PBEST ) # A 1 8 1 PBEST )*•••#  A (PBEST I PBEST ) I OUTPUT 
C PROD  a PRODUCT (I-A (PI PBEST )««8)  FOR  PfI  TO  PBEST)  OUTPUT 
t RHO(I) # • • • iRHO(PMAX)  ■ NORMALIZED  CORRELATION  COEFFICIENTS!  OUTPUT 
C DIMENSION  X(.0#Y(N)#A(PMAX)#HHO(PMAX)  is  reouired  in  main  prooram 
INTEBER  PMAX#Pat5T»P 
DOUBLE  PRECISION  SA#SB 
DIMENSION  X(i)#Y(l)#AUI»RHO(i) 
lF(PMAXaST#a#aS8RT (N) ) PRINT  8#  PMAX#N 
g Format (/t  pmax  a*#u#*  is  too  laroe  for  number  of  data  points  n a* 

B#  IB) 

C COMPUTE  MEAN 

Slao. 

uo  I lai.N 
I Siasi+X(l) 

Sissi/N 

C SUBTRACT  mean#  And  SCALE  TO  UNIT  VARIANCE 
S2a0. 

DO  4 Iai,N 
X(XJ8X(I)-SI 
3 SiaS24X<I)**8 

S28SQRT<sa/(N-l,n 
Tsl,/S2 
DO  5 1*1, N 
A ( I ) »X ( I )*T 
b Y(I)ax(I) 

C BCQiU  RECURSION 
P*U 

PHOOUCal, 

AICMINao, 

PbESTaO 
PHUOal . 

Psp+I 

CALCULATE  CR0SS-8AIM  £<J.  IBS 
SA=0. DO 
SttCO.UO 
LsP+1 

Du  7 I*L#N 
sabsa+x  t i )*y  ( i**i) 

SdaSD^X ( I ) **i* Y ( 1-1 ) **2 
tia2,*8A/SB 

hruducsproduc*  ( i , -i#*e ) 

CALCULATE  FILTER  COEFFICIENTS)  EOS.  140414A.  STORE  IN  

XiH)BO 

IKlP.EO.il  80  TO  b 

L=P/2 

Du  V IsitL 
TsX ( I )“G*X(P“I ) 
xiP-i)sx(P-i)-e*x(i) 

A(I)»T 

CALCULATE  NORNALUED  CORRELATION  COEFFICIENT  I EO.  149 
TsA(P) 

IFlP.Eu.l)  00  TO  14 
LsP-1 

Uu  lb  lsi#L 
IE  7sUACI)*RhO(P-I) 

14  PhU<P)*"i 

J-4 


■>  v . A .*#■■  • ....  ■ 
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C CALCULATE  AKAIAE'S  INFORMATION  CRITERION!  ESS. 
RtLEKRs  1 1 .-&*• ) *$N4L  < Sb  ) / (B  . * <N-P>  i 
A XCsLOO ( RELERR) ♦* . *FLOAT ( P ) / ( N-P ) 
lp( AlC.Ot. AXCMlNI  SO  TO  10 
AICMINXAXC 
HtttSTaP 
Pk00»PR0i)UC 
Du  11  lal»P 
11  Au)aX(I) 

10  1F(P.ES«PMAX>  SO  TO  IS 

C uPQATfc  FORoaRO  ANO  BACKWARD  SEQUENCES l EQ.  IBS 
LaP+1 

UO  IB  IaN*S-*“l 
YaX(l)-6*YtI-l) 

IB  XU)»T 


Y(P)a« 

00  TO  A 

16  Y(PMAX)aQ 

IF (PBEST rffcQ.PMAX!  00  TO  4 
C COMPUTE  EXTRAPOLATEO  NORMALIZED  CORRELATION 
C COEFFICIENTS  FROM  PBEST*!  TO  PMAXI  la.  105 
LaPBEST+1 
00  IT  PbL»PMAX 
A(P)>0. 

TaO, 

00  10  ui . PBEST 

16  TaWA(I)«RHO(P*l> 

17  NHO(P)aT 

4 X(N)mSl 

Y ( N ) aSB 
RETURN 

end 


166620* 


j 

'1 


3 


f. 

6 


[ 

I 


SUBROUTINE  POWERS (PBEST »A»PROD» J*XX*Yt.CO) 

C THIS  SUBROUTINE  COMPUTES  THE  FRACTIONAL  POWERS  IN  BANDS  1/<J*0ELTA)I  CO.  J-4 
C PBEST  a BEST  ORDER  OF  FILTER!  XNTEOER  INPU* 

C A ( 1 )».,.* A ( PBEST ) a FILTER  COEFFICIENT  ARRAY I INPUT 
C PROD  a PRODUCT (l-A(P) *»8>  FOR  Pal  TO  PBEST*  INPUT 
C J a SUE  OF  FFT  ( J/BUsNUmBEr  OF  FREQUENCY  POINTS)!  INTEQER  INPUT 
C XX  a AUXILIARY  ARRAY  ON  INPUT 
C XX(1)',,,'XX(U/B+1)  a FRACTIONAL  POWERS  ON  OUTPUT 
C YY  a AUXILIARY  ARRAY!  SCRATCH  INPUT 

C C0(1)',,"C0(JS4+1)  a QUARTER  COSINE  TABLE  FOR  FFT)  INPUT 
C DIMENSION  XX(U)’YY(U)#C0(U/4*1)  IS  REQUIRED  IN  MAIN  PROQRaM 
C DIMENSION  AIPMAX)  IS  REQUIREO  IN  MAIN  PROQRAM i WHERE  PMAX.QE. PBEST 
integer  pblst 

DIMENSION  A(l) tXX(l) rYY(l) tCO(l) 

FaPR004B«/J 

XX(l)al, 

YVdlao. 

DO  1 lai, PBEST 
XX(l+l)3»A(X> 

1 YY(I*l)aO. 

LaPBEST+Q 
OU  a UL.U 
XX(I)aO, 
k YY(I)aO, 

Lal,44B7*L0Q(J)*.8  0 L0#2(J) 

CALL  M«lFFT(XXiYY.CO»L#-l> 

Lbu/2+1 
DO  3 I*l»L 

3 XX(I)aF/|XX(I)**B+YY(I)MB) 

RETURN 

END 


* 
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1NFU1  DATA) 


•,b3*B?u96*0i 

.113*0011*04 

.18*08*18 *08 

,7+9436*6+01 

-,227871+0+01 

-.81771*14*01 

.67*0*6*0*01 

.80610847*08 

•1*86*686+08 

,1060377**01 

-,JbB*ou«,3»0* 

•<*11**648*04 

-.30**3711*08 

.289*1067*01 

.3676*337*00 

.84+01083+04 

.3607**33*04 

•■78**80*8*01 

-.*13*6869+08 

-.619*3+76*09 

-,J*++b#7**04 

.37*11677*01 

.*667120**08 

* 302*01*9+04 

.3*106104*02 

-,4+300to7*0l 

•.3*3*1616*08 

-.62051404*02 

•.37+37309+04 

-.68*96+64+01 

.I5o8i4t9+04 

.*3686*48*04 

.*1103811*04 

.437*6403+02 

36*1+612*01 

-.8637143**04 

•i*036«bl3*04 

-.3838**68*08 

-.1817663**04 

.136+8936+04 

,06701371*04 

.*73*3318*04 

.3376*6*8*04 

.16*7*633+01 

-.3*617383+12 

•.8**73361*04 

-.00417300*04 

-.18068873*08 

.423**9+1+04 

.+667670**02 

.+6*66046+04 

.410*7*18*08 

-.11764*31*01 

-.3201*0*6+02 

-.32151017+04 

-.l7»8*UkO*04 

.1*807*60*01 

.1766*800*02 

,46891872+02 

.29831*72+02 

,4k+5V073+0k 

.07883*77*01 

-.1*018278*01 

-.32*36678*02 

37616888*04 

-.8+*3..kbO*04 

-.*3786468*0(1 

.43166670*04 

.38*15088+04 

.316+5704+04 

.147*1387*04 

-.1106**74*08 

-.461*87***01 

-.31179**7+04 

-.19881729*08 

-.1*1060+5+01 

• 14J0t»63**04 

.18***507*08 

.i*Sl*t6++0» 

.13328086+02 

.33383348*00 

-.14632438*04 

-.38074808+00 

-.1812*2*3+04 

-.11167069+02 

.13083073*04 

■*u*l0*0**02 

.3*112601+04 

.++7700*3+08 

-.376+1+42*01 

-,5bb06b50+04 

-.30*87314*08 

-.37683633+04 

-.118117+7+00 

,89441863+02 

.65*0374**04 

.*33603***04 

.60*09774+01 

-.32971*00+04 

-.793+8*89*02 

htAN  & .11733613*00 

xi+i.uahu  ut. vim r io,~i  x .i;36too7+oa 


I'ulbT  .1  + 


hi+oicrm  rikTth  eufcFkicimi >i 

•47038131*01  -.3784Jblb*Cl 

.46*06147+01 

-.+3677754*00 

^ko(iOCT(l-*ll'l.*4)  1 , 

22986847-83 

NUh+iAt  1 41U  CUHKEL+TIUN 
.73304340+00 
-.43187107*80 

C08F91CIIN7SI 
.2110+268-01 
. 1*061 780+00 

-.S*61«+e*+('0 

,76830693+00 

-.*3601*70*00 

.73bJ+707+u0 

-.777ull.27*00 

,3or+362*+lin 

BAD  DATA  POINTS  '(SUBSECTION  5. 2) 


C SPECTRAL  ESTIMATION  FOR  BAD  DATA  POINTS  USER I CHANSE 

C LINE  17  ANO  REPLACE  LINES  |2-3t  AND  41-Hi 
C N i NUMBER  OF  DATA  POINTS 
C XUImioKN)  b INPUT  DATA 
C UMAX  v MAXIMUM  NUMBER  OF  BAD  DATA  POINTS 

c a ■ ACTUAL  NUMBER  QP  BAD  DATA  POINTS  (MUST  HAVE  B.LE.BMAX) 

C M(1)»,,,»M(B>  b LOCATIONS  OF  BAD  DATA  POINTS 
c pmax  » maximum  order  of  filter 
C PBEST  a BEST  ORDJCR  OF  FILTER 

C A(l) r • • •» A (PBEST)  a PREDICTIVE  FILTER  COEFFICIENTS 
C PROD  ■ PRODUCT U-A(P>**B)  FOR  Pal  TO  PBEST 

C RHO(l) > . • • iRHO(PmAX)  b normalized  correlation  coefficients 

C U b SIZE  OF  FFT  (MUST  BE  A POKER  OF  i) 

C XX 1 1) » , . • t XX( J/2+L)  a FRACTIONAL  POBERSt  FROM  DC  TO  NTOUIST  FREBuENCV 
C C0(l).,,,fC0(U/4+t)  B ouap.ter  cosine  table 
C Y*  YY*  AND  IP  are  REfiUlRED  AUXILIARY  ARRAYS 

parameter  nb  ioo»  bmaxb  a,  pmaxbxq,  0«ma»  unajssn 

1NTE0ER  9* PBEST 

DIMENSION  X(N)  #Y(N)  t A(PmAX)  «RH0(PMAX)  »XX(+I)  »YY  (J)  »C0(J41i 
DIMENSION  M (BMAX) r IP (N) 

C INPUT  DATA  IN  XU)m..*X(N> 

DEFINE  IRANDBI *8**  18+  ( ( 1-SX9N  ( 1 , 1 *0***6 > ) /* ) *01*369736 36,’ 

DEF INC  RaNDbFlO AT ( l ) /343»«736367 * 

136261 

NSTaHTsN+600  « MILL  DISCARD  INITIAL  HOO  POINTS 
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hi 


xxinao, 

xxiaiso, 

Xx(3)a0, 

XXI4>aO, 

Uu  11  LaSiNSTANT 
1*  IRANI) 

XX(L)82.7607*XA(L-l)-3,8l0S*XXlL-ai*8,653b*XXlL-3)« 

*0,V836«XaIL-4)+TINORM(RaND,811) 

CONTINUE 

00  18  XsliN 

XU)»XX(I+NSTART-N) 

PRINT  1» 

fokmat i/t  input  dataim 

PkINT  4,  ( X 1 1 ) » laliN) 

ENTER  Bt  AND  ENTEN  BAD  OATA  LOCATIONS  IN  MU)'»..'M(B) 
bob 
K(l)s3 
Mi2)a7 
M(3)sli 
M4)al2 
M(a)al9 

EVALUATE  PREDICTIVE  filter  coefficients 

call  BURSBO (N , PMAX » X * tt i M# IP » V » PBEST . A , PROD » RHO ) 
PRINT  9»  X (N) 

FORMAT  1/ 1 MEAN  b'*E14*6> 

PRINT  10i  T(N) 

FORMAT ( ' STANDARO  DEVIATION  s»,El3.Bl 

PRINT  2.  PBEST 

FORMAT (/!  PBEST  S*«X3» 

PRINT  3. 

FORMAT </!  PREDICTIVE  filter  COEFFICIENTS!') 

PRINT  4,  lA(I)«lal»pBEST> 

F0RMAT<SE80.B) 
print  b,  prod 

FghMAT (/!  PRODUCT (1-A<P)**8)  ='»E13.8) 

PRINT  t>» 

FORMAT!/'  NORMALISED  CORRELATION  COEFFICIENTS! • > 
PRINT  4,  (KhOU)*I«1#PMaX) 

CALL  QTRCOS(COiU) 

evaluate  fractional  powers 

call  POWtRS<PaeST»A,PROD*U»XXfYY»CO> 

Print  7, 

FORMAT!/!  FRACTIONAL  P0wER5» ') 

Lsii/2+1 

PRINT  B i (XX! I ) • I*1»L) 

FORMAT (2x»10El3,b) 

End 


subroutine  BUROBO ( N » PMAX  * X»B#M»IP»Y*  PBEST i A > PROD  * RHO ) it  2 FEB  1976 
THIS  SUBROUTINE  COMPUTES  ThC  PREDICTIVE  KILTER  COEFFICIENTS  FOR  B BAD  POINTS 
N a NUMBER  OF  DATA  POINTS!  XNTESER  INPUT 
PMAX  m MAXIMUM  ORDER  OF  FILTER I XNTESER  INPUT 
Xll)'X(2)'.,.#X(N>  l OATA  ARRAY  ON  INPUT!  ALTEREC  ON  OUTPUT 

ON  OUTPUT » X(l)’X(2)iit» tX(PMAX)  « A(l IPMAX) i A A<PMAX(PMAX) 

B a NUMBER  OF  BAD  OATA  POINTS!  XNTESER  INPUT 

M(l) »M(8) • ,t% «M(B)  B LOCATIONS  OP  BAD  DATA  POINTS!  XNTESER  INPUTS 
THESE  LOCATIONS  MUST  BE  DISTINCT  AND  LIE  IN  THE  RAN9E  Cl'Nl 
IP(l)iXP(2)i...*XP(N)  a AUXILIARY  ARRAY!  SCRATCH  INPUT 
Y(l) i Y{8) » •/•! Y(N)  *>  AUXILIARY  ARRAY!  SCRATCH  INPUT 

ON  OUTPUTt  Y!1)'Y(2)m.«»Y  (PMAX)  a Mill) . A (PMAX  IPMAX) 

ON  OUTPUT  * X(N)  a MEAN*  AND  YIN)  a STANDARD  DEVIATION  QP  INPUT  DATA 
PBEST  a BEST  ORDER  OP  FILTER!  XNTESER  OUTPUT 
All) » A(8) t , ,t r AlPBEST)  a PREDICTIVE  FILTER  COEFFICIENT  ARRaY  a 
A U I PBEST ) *8(81 PBEST )#..., A (PBEST I PBEST ) f OUTPUT 
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I 

i 


f 

f 


i. 

i 

f 

r 


C PrtOO  c PRODUCT ( 1-A (PI PBEST ) **a ) POP  Pal  TO  PBESTI  OUTPUT 

0 RHO(l) • «•« #RHO{PMAX)  ■ NORMALIZED  CORRELATION  COEPPICIENTSl  OUTPUT 

C DIMENSION  X(N)*Y(N)»A<PMAX)*RHO<PMAX)  IS  RESUMED  IN  MAIN  PROGRAM- 
C DIMENSION  M(BMAX)iIPIN)  IS  REflUIRCD  IN  MAIN  PROORAM 
INTEOER  PMAX * B r PBEST iPidP 
DOUBLE  PRECISION  SA.SB 
DIMENSION  X(l).MU)  ilP(l)  *V{l)»A(lliRftO(l) 
iP(e.OT.o)  oo  to  ai 

CALL  BURO(N>PMAXiXfV«PbESTiAiPRODiRHO) 

RETURN 
XI  LSN-B 

IF(PMAX.OT.S.*SBRT(L>>  PRINT  2r  PMAX.L 
k FORMAT!/)  PMAx  a'.l4,'  IS  TOO  LARBE  FDR  NUMBER  OF  800D  DATA  POINTS 

S N»B  » tSi 

C SET  UP  IP  ARRAY  FOR  P«OI  E«.  ITS 
oo  22  ibi.n 
22  iPtnai 

OU  as  Lsl.B 
laM(L) 

kS  IP(l)aO 

C COMPUTE  MEAN  OF  OOOD  UATA  POINTS 
S1B0. 

00  1 lBltN 

IP(IP(I),E8,0)  00  TO  1 
SlaSl+X(l) 

1 continue 

SlaSI/INwBI 

C SUBTRACT  MEANi  and  scale  to  unit  variance*  fur  bood  data  points 
SitBO. 

DO  S I«1»N 

IFIIP(I) ,EQaO)  bo  TO  3 
X(I)*X(IJ-S1 
SkB6a+X(i)**a 
3 CONTINUE 

sxbsqrt ( sa/ ( N-B-l • ) ) 

Tal,/sa 
DO  0 Iai,N 

IF( IP ( I ) ,E8,0)  SO  TO  8 
X(!)aX(!)«T 
YU)«X(lj 
b CONTINUE 

C BE8IN  RECURSION 
PaO 

PHODUCsi* 

AICMINao, 

PdESTsO 

PROual. 

b Pzp+1 

C UPUArt  IP  ARRAY)  £0*  173 
DU  24  Laird 
I=M(L>*P 

1FU.BT.N)  60  TO  24 
IP(I)aO 

*4  CONTINUE 

BPso 
LSP+1 

DO  28  1sl«N 
as  BPS8P41-JP!D 
KsN-P-bP 

iFix.LT.ab)  print  tt>  k,p 

*6  F0RmAT(/1  number  OF  VALID  ERROR  POINTS  IS  0NLY'rI3»*  FOR  ps«H3) 

C CALCULATE  CROSS-OAINI  EQ.  1*3 
SABO, 00 


j-e 


r 

l 

i 


LL. 


>■  iuudi  ^ - -i.  m . i L . . JiA.  . ^ 


/ 
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i 


» 


7 

C 


6 

C 

tt 


18 

14 

C 


11 

10 

c 


12 


16 

L 

t 


16 

17 

4 


*•*0,00 

L*P+1 

00  7 S*liN 

iFUPtx>,Eu,o>  «o  to  7 

*A*SA*Xu>*V'I~1> 

Sb*8B+X  ( l ) 4*8+7  ( 1-1 ) +*2 

CONTINUE 

0*2, *S A/SB 

PROOUC*PROOUC  * ( 1 , -0*6 ) 

CALCULATE  FILTER  COEFFICIENTS!  EOS,  1S80&96.  STORE  IN  XU) 
XiP)*B 

IF(p,e«,l)  «C  TO  0 
LsF/2 

00  9 1*1,  L 

TSAf l)-6*X{p-i) 

x<p-n*xiP-i)-o+xu) 

XU)*T 

CALCULATE  NORMALIZED  CORRELATION  COEFFICIENT!  EO,  149 
T*A(P) 

IFCH.Ew, 1)  00  TO  14 
L*P-1 

DO  IQ  XalfL 
T*T+XU)+RHO<P-I) 

RHO(P)*T 

CALCULATE  AKAXKE'S  INFORMATION  CRITERlONf  EOS,  1940202 
HfcLERHs U .-0*0 ) +SN8L ( SB ) / ( 2 , +K ) 

A 1C*L00 ( hELERR ) +2 , *FLO AT (P ) /K 

IFUIC.OE.AICMIN)  00  TO  10 

AICMIN«A1C 

PdESTaP 

PHOO*PHOOUC 

00  ll  Isl»p 

AU»*XU) 

XF(P,S6,PMAX>  00  TO  16 

UPDATE  FORWARD  AnD  BACKWARD  SEQUENCES*  E0.101 
LsP+1 

DO  12  IeN»L»-1 

IF  IIP(I) ,EO,0)  00  TO  12 

T=A(I)-0*V(I-i) 

7(l)*Y(l«i)-0*X(I) 

XU)*T 

continue 

7 (P)sO 
60  TO  6 
ViPMAX)sg 

IF (puEST ,tu,PMAX)  60  TO  4 
CONFUTE  EXTRAPOLATED  NORMALIZED  CORRELATION 
COEFFICIENTS  FROM  PBEST+1  TO  PHAXI  E9.  168 
LsPBEST+1 
Do  17  PsL*PNAX 
A IP ) *0, 

7*0, 

DO  16  I*1»PbEST 
TsT+Atl)*RHO(P-I) 

Rhu(P)*T 
X ' N)*S1 
Y (n)*S2 

return 

end 


X(P) 
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SUBROUTINE  POWERS (PBEST » A . PROOiJtXX .YY.CO) 

C THIS  SUBROUTINE  COMPUTES  THE  FRACTIONAL  POWERS  IN  BANOS  1/(U*DELTA)I  KB.  7-R 

C PBEST  s BEST  ORDER  OF  FILTER)  INTEBER  INPUT 

C A(l)  t .. ..A(PBCST)  » FILTER  COEFFICIENT  ARRAY  I INPUT 

C PROD  B PRODUCT U-A IP) *•*)  FOR  PB1  TO  PBEST I INPUT 

CU«  SIZE  OF  FFT  (U/B*I«NUMBER  OF  FREQUENCY  POINTS)!  INTEBER  INPUT 

C XX  ■ AUXILIARY  ARRAY  ON  INPUT 

C XXU)',t.’XX(U/B«l>  > FRACTIONAL  POWERS  ON  OUTPUT 

C YY  ■ AUXILIARY  ARRAY!  SCRATCH  INPUT 

C COa).M.’CO(U/*«l)  > QUARTER  COSINE  TABLE  FOR  FFT!  INPUT 
C DIMENSION  XX(U)'VV(U>tC0(J/4+l)  IS  REQUIRED  IN  MAIN  PROORAM 
C DIMENSION  A(PMAX)  IS  REQUIRED  IN  MAIN  PRQBRAMt  WHERE  PMAX.BE, PBEST 
INTEBER  PBEST 

DIMENSION  A(l) rXX(l) t YY(1) tCO(l) 

FBpROO«ti/U 

XX<i)«i, 

TY(1)bO, 

DO  t I Bit PBEST 
XX(ln)l-A(X) 

1 VY(I>1)bO. 

LQPBE8T4* 

DO  8 X«LtJ 

xxmso, 
a YY«n*o, 

L> 1 . 44B7PL0Q ( J ) * • 8 LOOa(J) 

CALL  MKLFFT ( XX  t Y Y i CO  t L » »1 ) 

LBJ/B41 
00  3 IiltL 

3 XX(I)aP/(XX<I)**a4YY(I)B*a) 

RETURN 

ENO 
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